Method for laying forest fixed large sample plot by using RTK (Real-Time Kinematic)

Abstract

The invention discloses a method for arranging a forest fixed large sample plot by applying RTK, relates to the technical field of forestry investigation, and aims to solve the technical problems that manual point arrangement and tree species identification mark point arrangement are adopted for current forest fixed large sample plot point arrangement, time and labor are wasted, and efficiency and precision cannot be effectively improved. An unmanned aerial vehicle is used for carrying a laser radar, airborne laser scanning is carried out above an original terrain, high-precision position coordinates are obtained through GNSS / INS combined positioning, and point cloud data are collected; s2, establishing a three-dimensional digital model; processing the point cloud data, filtering noisy points, removing isolated points, generating a digital elevation model based on the processed point cloud data, and generating a forest fixed large sample plot center point in the digital elevation model; s3, field sample plot planning; and dividing the layout area of the fixed large sample plot into a plurality of field sample plots. The method has the advantage that the positioning precision and the sample plot laying efficiency are improved.

Description

Technical Field

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

[0002] With the deepening of research on global forest ecosystems, 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 forest ecosystem research. Dynamic monitoring of forest resources mainly includes real-time and regular monitoring of resource parameters such as forest stock, area, distribution of tree species and canopy.

[0003] Traditional forest ecology plot research is to establish ground plots and apply corresponding geostatistical methods to conduct small-scale forest resource surveys to obtain characteristic data for each type of tree in the plot, and then establish a model to simulate the dynamic changes of forest resources. The layout method of fixed forest plots is mainly based on traditional forest resource survey methods. This method requires a lot of manpower and time costs, and is not easy to carry out in large ecosystems.

[0004] The use of remote sensing image analysis methods to dynamically monitor forest resources can save manpower and time costs of forest surveys and improve efficiency. It is currently the main method. However, in remote sensing analysis of forest resources, it is difficult to obtain accurate ground verification data using traditional methods, and there are different opinions on methods such as image extraction. Especially when dealing with complex terrain, such as hilly terrain, the analysis of remote sensing images cannot effectively obtain the true value of the image. Therefore, it is very necessary to use the research of forest fixed ecological plots to strengthen the dynamic monitoring of forest resources.

[0005] At present, there are two main ways to set up forest fixed ecological plots: manual point setting and tree species identification mark setting. Manual point setting requires reference to previous data and actual terrain conditions, representativeness of sample points, rich experience and professional knowledge, and consumes a lot of manpower and time; tree species identification mark setting requires manual visual markers, which need to be particularly obvious objects that can be identified from the air, and require certain area, shape and other characteristics, requiring special investment of manpower and time, and repeated field identification marking, which is time-consuming and labor-intensive, and cannot effectively improve efficiency and accuracy. In view of this, we propose a method for setting up large forest fixed plots using RTK. Summary of the invention

[0006] The purpose of the present invention is to provide a method for laying out large fixed forest plots using RTK, so as to solve the technical problem that the current large fixed forest plots are laid out manually and with tree species identification marks, which is time-consuming and labor-intensive and cannot effectively improve efficiency and accuracy.

[0007] In order to solve the above technical problems, the present invention provides the following technical solutions: a method for deploying a large forest fixed sample plot using RTK, comprising the following steps:

[0008] S1, point cloud data collection;

[0009] Use a drone equipped with a laser radar to perform airborne laser scanning above the original terrain, obtain high-precision position coordinates through GNSS / INS combined positioning, and collect point cloud data;

[0010] S2, establishment of three-dimensional digital model;

[0011] Processing the point cloud data, filtering out noise points, removing isolated points, and generating a digital elevation model based on the processed point cloud data, and generating a center point of a large forest fixed sample plot in the digital elevation model;

[0012] S3, Field plot planning;

[0013] Divide the layout area of ​​the fixed large sample plot into several field sample plots, and determine the scope of the field sample plot and the coordinates of the center of the field sample plot;

[0014] S4, obtaining the center of the drone's outward expansion route;

[0015] Based on the center coordinates of the actual sample plot, the center coordinates of the drone's outward expansion route are calculated, and the outward expansion flight altitude is determined according to the area of ​​the actual sample plot and the height difference between the highest point and the lowest point in the fixed large sample plot to generate the outward expansion flight route;

[0016] S5, obtaining the coordinates and elevations of the expanded sample points;

[0017] Use drones to perform airborne laser scanning along the outward expansion flight route to obtain the coordinates and elevations of the outward expansion sample points;

[0018] S6 lays out fixed large sample plots;

[0019] According to the coordinates and elevations of the expanded sample points, the coordinates of the control sample points of the fixed large sample plot are calculated in reverse, and the fixed large sample plot is arranged according to the coordinates of the control sample points.

[0020] Preferably, the method for filtering out noise points from point cloud data is based on a statistical filtering algorithm, and the algorithm formula is as follows: Suppose a point P in the point cloud data i (x i ,y i ,z i ), calculate the average distance to the surrounding k neighboring points And the average distance of all points is μ d and standard deviation σ d ;

[0021]

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

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

[0024] Preferably, in S1, the unmanned aerial vehicle is a fixed - wing unmanned aerial vehicle, and the digital elevation model generation method is based on the construction of an irregular triangular 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 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:

[0025]

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

[0027] Map the digital ground points to the plane of the digital ground elevation model to establish a digital ground model. Thinning the original terrain in the three - dimensional digital model within the area, and adding the digital height of the digital ground model to the digital ground points in the thinned digital ground model.

[0028] Preferably, the method for planning field sample plots for the three - dimensional digital model of the original terrain and calculating the center of the UAV's outward - expanding flight path by reversing from the center of the field sample plot includes:

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

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

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

[0032] a4. Based on the center coordinates of the field sample plot, obtain the center coordinates of the UAV's outward expansion route.

[0033] Preferably, the method of back-calculating the center of the drone's outward expansion route from the center of the actual sample plot is to back-calculate the center of the fixed large sample plot from the coordinates of the center of the actual sample plot within the actual sample plot range, and assume that the coordinates of the center of the actual sample plot are (X si ,Y si ), where i = 1, 2, ..., m, m is the number of field plots. First, the average value of the center coordinates of all field plots in the X direction and the average value in the Y direction are calculated. The formula is as follows:

[0034]

[0035] Then, As the initial estimate, based on the coordinates of the expanded sample points (X spj ,Y spj ,H spj ), j = 1, 2, ..., n, n is the number of external expansion sample points, and weighted least squares adjustment optimization is performed;

[0036] Assume that the distance from the expanded sample point to the initial estimated center is Construct the error equation:

[0037]

[0038] where v Xj 、v Yj is the observation residual in the X and Y directions, δX and δY are the correction values ​​of the center coordinates in the X and Y directions, Δd j It is the correction number of the distance from the expanded sample point to the optimized center coordinate;

[0039] According to the least squares principle Construct the normal equation to solve δX and δY, where p j The weight of the external sample point is determined according to the accuracy of the point, and the normal equation is constructed to solve δX and δY;

[0040] Finally, the center coordinates of the large sample site (X c ,Y c ) is calculated as:

[0041] The center coordinates of the fixed large sample plot are superimposed on the three-dimensional digital model of the original terrain to obtain the center coordinates of the outward expansion route;

[0042] Assume that the center coordinates of the actual sample plot are (Xsi ,Y si ), the offset of the center of the fixed large plot relative to the center of the actual plot in the X direction is ΔX, and the offset in the Y direction is ΔY, then the coordinates of the center of the UAV's outward expansion route (X o ,Y o ) is calculated as:

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

[0044] Preferably, the method for determining the outward expansion flight altitude of the drone through the outward expansion route center 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 outward expansion route center according to the height difference, and set the flight altitude to the elevation value from the lowest point in the fixed large sample plot to the outward expansion route center according to the distance;

[0045] Let the coordinates of the center of the actual plot be , and let the elevation of the highest point in the fixed large plot be H max , the lowest point elevation is H min , the horizontal distance from the lowest point in the plot to the center of the outward expansion route is D, then the outward expansion flight altitude H filght The calculation formula is: Where (ΔH) = H max -H min .

[0046] Preferably, the method for calculating the coordinates and elevations of the expanded sample points through the expanded sample points includes:

[0047] b1. Using the center coordinates of the fixed large sample plot to overlay the three-dimensional digital model of the original terrain, aligning the center coordinates of the fixed large sample plot with the outward expansion flight route to obtain the outward expansion flight route;

[0048] b2. Use a drone to fly at the altitude value according to the outward expansion flight route to perform outward expansion scanning and obtain the coordinates and elevation of the outward expansion sample points;

[0049] Assume that the center coordinates of the extended flight route are (X o ,Y o ), the offset distance of the expanded sample point relative to the center of the expanded flight route in the X direction is ΔX spj , the offset distance in the Y direction is ΔY spj , the radius of the earth is R, the angle between the line connecting the outward expansion point and the center of the outward expansion flight route and the east direction is θ, then the coordinates of the outward expansion point (X spj ,Y spj ) is calculated as:

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

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

[0052] Where ΔX spj and ΔY spj According to the set expansion distance interval and the expansion point sequence number, let the expansion step length be d, and the sequence number of the nth expansion point in this direction be n, then ΔX spj =n×d×cosα,ΔY spj =n×d×sinα;

[0053] Assume that the flight altitude is H flight , the center elevation of the extended flight route is H o The vertical height deviation of the expanded sample point relative to the center of the expanded flight route is ΔH spi , then the elevation of the extended sample point H spi The calculation formula is: spi =H o +H flight +ΔH spi .

[0054] Preferably, the outward expansion flight route is outward expansion from the center of the outward expansion flight route, and the scanning points are the outward expansion sample points. The method for calculating the coordinates and elevations of the outward expansion sample points through the outward expansion sample points is to obtain distance deviations based on the outward expansion sample points and the center of the outward expansion flight route, determine the heading distances of adjacent outward expansion sample points based on the distance deviations, and determine the longitude and latitude of the outward expansion sample points according to the heading distances;

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

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

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

[0058] Preferably, 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:

[0059] Assume that the coordinate set of the external expansion sample points is {(X spi ,Y spi ,H spi )} i n ＝1 , 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, and the fitting parameters are calculated using the following formula:

[0060]

[0061] Solve for a, b, c;

[0062] Assume that the fixed plot is a rectangle with a length of L and a width of W, and fit the coordinates of the center point of the plot on the plane (X c ,Y c ,H c ) is used as the reference, and m control sample points are set at equal intervals in the long direction i. In the width direction j, k control sample points are set at equal intervals. Then the coordinates of the control sample point in the jth row and the ith column (X cpij ,Y cpij ,H cpij ) is calculated as follows:

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

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

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

[0066] The pole markers are set at the control sample points by using the GNSS instrument to obtain the coordinates of the pole markers, and the layout direction of the pole markers is set to be consistent with the fixed sample plot.

[0067] Preferably, the coordinate accuracy calculation of the pole marker layout is also included, and the accuracy values ​​of the pole marker coordinates and the control sample point coordinates are calculated through the external expansion sample points and the external expansion flight route. The calculation method is: let the theoretical coordinates of the control sample point calculated through the external expansion sample point be (X ct ,Y ct ,H ct ), the actual pole coordinate is (X rb ,Y rb ,H rb ), the precision threshold is ε, then the precision value P is calculated as:

[0068]

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

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

[0071] 1. The present invention uses a laser radar mounted on an unmanned aerial vehicle to collect point cloud data, and then reversely calculates the center coordinates of the unmanned aerial vehicle's outward expansion route according to the center coordinates of the actual sample plot, determines the outward expansion flight altitude, obtains the coordinates and elevation of the outward expansion sample points, and finally reversely calculates the coordinates of the control sample points of the fixed large sample plot, which greatly reduces manpower input and is beneficial to improving positioning accuracy and sample plot layout efficiency.

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

[0073] 3. The present invention can also obtain rich terrain and land object information through point cloud data collected by drones. Compared with the limited field measurement data or remote sensing image data that are difficult to accurately interpret that may be relied upon in the prior art, point cloud data can more comprehensively reflect the actual situation of the sample site and provide more complete data support for forest resource surveys and research. DETAILED DESCRIPTION

[0074] In order to facilitate those skilled in the art to understand the technical solution of the present invention, the technical solution of the present invention is now further described.

[0075] A method for deploying a large forest fixed plot using RTK, comprising the following steps:

[0076] S1, point cloud data collection;

[0077] Use a drone equipped with a laser radar to perform airborne laser scanning above the original terrain, obtain high-precision position coordinates through GNSS / INS combined positioning, and collect point cloud data;

[0078] S2, establishment of three-dimensional digital model;

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

[0080] S3, Field plot planning;

[0081] Divide the layout area of ​​the fixed large sample plot into several field sample plots, and determine the scope of the field sample plot and the coordinates of the center of the field sample plot;

[0082] S4, obtaining the center of the drone's outward expansion route;

[0083] Based on the center coordinates of the actual sample plot, the center coordinates of the drone's outward expansion route are calculated, and the outward expansion flight altitude is determined according to the area of ​​the actual sample plot and the height difference between the highest point and the lowest point in the fixed large sample plot to generate the outward expansion flight route;

[0084] S5, obtaining the coordinates and elevations of the expanded sample points;

[0085] Use drones to perform airborne laser scanning along the outward expansion flight route to obtain the coordinates and elevations of the outward expansion sample points;

[0086] S6 lays out fixed large sample plots;

[0087] According to the coordinates and elevations of the expanded sample points, the coordinates of the control sample points of the fixed large sample plot are calculated in reverse, and the fixed large sample plot is arranged according to the coordinates of the control sample points.

[0088] In an embodiment of the present invention, the method for filtering out noise points in point cloud data is based on a statistical filtering algorithm, and the algorithm formula is as follows: Suppose a point P in the point cloud data i (x i ,y i ,z i ), calculate the average distance d to the surrounding k neighboring points i , and the average distance of all points μ d and standard deviation σ d ;

[0089]

[0090] like Then the decision point P i The noise points are filtered out, where n is the set threshold coefficient, which is 2-3;

[0091] The method of removing isolated points from point cloud data is based on cluster analysis algorithm. For any point P in the point cloud data, i , calculate the number of neighboring points N(P i ,r), if N(Pi , if r < m, then mark point P i as an isolated point and remove it, where m is the minimum clustering point number threshold.

[0092] In an embodiment of the present invention, in S1, the unmanned aerial vehicle is a fixed-wing unmanned aerial vehicle, and the digital elevation model generation method is based on the construction of a triangulated irregular network (TIN), and the generation method is as follows: First, perform Delaunay triangulation on the processed point cloud data to construct a triangulated irregular network. Let any triangle T abc in the triangulated 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:

[0093]

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

[0095] Set digital ground points and a digital ground elevation model plane 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.

[0096] The original terrain is a hilly terrain. 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.

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

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

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

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

[0101] a4. Based on the center coordinates of the actual sample plot, obtain the center coordinates of the UAV’s outward expansion route.

[0102] In an embodiment of the present invention, the area of ​​the field plot is 20m×20m. The selection of the side length of the field plot is based on a comprehensive consideration of the existing plot layout. A small plot area is prone to point data loss, a large plot area is prone to irregular plot shape, and a large field plot area affects the sampling accuracy in the later field survey. The size of the field plot area is based on the actual survey area and the terrain conditions of the actual survey area. The method of back-calculating the center of the drone's outward expansion route through the center of the field plot is based on the field plot range. The center of the fixed large plot is obtained by back-calculating the center coordinates of the field plot. The coordinates of the center of the field plot are (X si ,Y si ), where i = 1, 2, ..., m, where m is the number of field plots. First, the average value of the center coordinates of all field plots in the X direction is calculated. and the average value in the Y direction The formula is as follows:

[0103]

[0104] Then, As the initial estimate, based on the coordinates of the expanded sample points (X spj ,Y spj ,H spj ), j = 1, 2, ..., n, n is the number of external expansion sample points, and weighted least squares adjustment optimization is performed;

[0105] Assume that the distance from the expanded sample point to the initial estimated center is Construct the error equation:

[0106]

[0107] Where V Xj 、v Yj is the observation residual in the X and Y directions, δX and δY are the correction values ​​of the center coordinates in the X and Y directions, Δd j It is the correction number of the distance from the expanded sample point to the optimized center coordinate;

[0108] According to the least squares principle Construct the normal equation to solve δX and δY, where p j The weight of the external sample point is determined according to the accuracy of the point, and the normal equation is constructed to solve δX and δY;

[0109] Finally, the center coordinates of the large sample site (X c ,Y c ) is calculated as:

[0110] The center coordinates of the fixed large sample plot are superimposed on the three-dimensional digital model of the original terrain to obtain the center coordinates of the outward expansion route;

[0111] Assume that the center coordinates of the actual sample plot are (X si ,Y si ), the offset of the center of the fixed large plot relative to the center of the actual plot in the X direction is ΔX, and the offset in the Y direction is ΔY, then the coordinates of the center of the UAV's outward expansion route (X o ,Y o ) is calculated as:

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

[0113] In an embodiment of the present invention, the method for determining the outward expansion flight altitude of the drone through the outward expansion route center 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 outward expansion route center according to the height difference, and set the flight altitude to the elevation value from the lowest point in the fixed large sample plot to the outward expansion route center according to the distance;

[0114] Let the coordinates of the center of the actual plot be , and let the elevation of the highest point in the fixed large plot be H max , the lowest point elevation is H min , the horizontal distance from the lowest point in the plot to the center of the outward expansion route is D, then the outward expansion flight altitude H flight The calculation formula is: Where (ΔH) = H max -H min .

[0115] In an embodiment of the present invention, a method for calculating the coordinates and elevations of the expanded sample points through the expanded sample points includes:

[0116] b1. Use the fixed large sample site center coordinates to superimpose the three-dimensional digital model of the original terrain, align the fixed large sample site center coordinates with the external expansion flight route, and obtain the external expansion flight route;

[0117] b2. Use a drone to fly at the altitude value according to the outward expansion flight route to perform outward expansion scanning and obtain the coordinates and elevation of the outward expansion sample points;

[0118] Assume that the center coordinates of the extended flight route are (X o ,Y o ), the offset distance of the expanded sample point relative to the center of the expanded flight route in the X direction is ΔX spj , the offset distance in the Y direction is ΔY spj, the radius of the earth is R, the angle between the line connecting the outward expansion point and the center of the outward expansion flight route and the east direction is θ, then the coordinates of the outward expansion point (X spj ,Y spj ) is calculated as:

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

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

[0121] Where ΔX spj and ΔY spj According to the set expansion distance interval and the expansion point sequence number, let the expansion step length be d, and the sequence number of the nth expansion point in this direction be n, then ΔX spj =n×d×cosα,ΔY spj =n×d×sinα;

[0122] Assume that the flight altitude is H flight , the center elevation of the extended flight route is H o The vertical height deviation of the expanded sample point relative to the center of the expanded flight route is ΔH spi , then the elevation of the extended sample point H spi The calculation formula is: spi =H o +H flight +ΔH spi .

[0123] 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 height 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 through 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;

[0124] 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 dh , the radius of the earth is R, then the latitude and longitude of the expanded sample point are calculated as follows:

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

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

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

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

[0129]

[0130] Solve for a, b, c;

[0131] Assume that the fixed plot is a rectangle with a length of L and a width of W, and fit the coordinates of the center point of the plot on the plane (X c ,Y c ,H c ) is used as the reference, and m control sample points are set at equal intervals in the long direction i. In the width direction j, k control sample points are set at equal intervals. Then the coordinates of the control sample point in the jth row and the ith column (X cpij ,Y cpij ,H cpij ) is calculated as follows:

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

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

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

[0135] During the field deployment process, pole markers are set at the control sample points through GNSS instruments to obtain the pole marker coordinates. The number of pole markers set must be no less than 3 to be used for sample plot deployment. The deployment direction of the pole markers is set to be consistent with the fixed sample plot to ensure deployment efficiency. The number of pole markers set must be no less than 3 to facilitate determination of the deployment position of the fixed sample plot.

[0136] In an embodiment of the present invention, the coordinate accuracy calculation of the pole marker arrangement is also included. The accuracy values ​​of the coordinates of the pole marker and the coordinates of the control sample points are calculated by the external expansion sample points and the external expansion flight route. The calculation method is: let the theoretical coordinates of the control sample points calculated by the external expansion sample points be (X ct ,Y ct ,H ct ), the actual pole coordinate is (X rb ,Y rb ,H rb ), the precision threshold is ε, then the precision value P is calculated as:

[0137]

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

[0139] The embodiments of the present invention disclose preferred embodiments, but are not limited thereto. A person skilled 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 deviate from the spirit of the present invention, they are all within the protection scope of the present invention.

Claims

1. A method for deploying large forest fixed plots using RTK, characterized in that: The following steps are involved: S1, point cloud data collection; Use a drone equipped with a laser radar to perform airborne laser scanning above the original terrain, obtain high-precision position coordinates through GNSS / INS combined positioning, and collect point cloud data; S2, establishment of three-dimensional digital model; Processing the point cloud data, filtering out noise points, removing isolated points, and generating a digital elevation model based on the processed point cloud data, and generating a center point of a large forest fixed sample plot in the digital elevation model; S3, Field plot planning; Divide the layout area of ​​the fixed large sample plot into several field sample plots, and determine the scope of the field sample plot and the coordinates of the center of the field sample plot; S4, obtaining the center of the drone's outward expansion route; Based on the center coordinates of the actual sample plot, the center coordinates of the drone's outward expansion route are calculated, and the outward expansion flight altitude is determined according to the area of ​​the actual sample plot and the height difference between the highest point and the lowest point in the fixed large sample plot to generate the outward expansion flight route; S5, obtaining the coordinates and elevation of the expanded sample points; Use drones to perform airborne laser scanning along the outward expansion flight route to obtain the coordinates and elevations of the outward expansion sample points; S6. Lay out fixed large sample plots; According to the coordinates and elevations of the expanded sample points, the coordinates of the control sample points of the fixed large sample plot are calculated in reverse, and the fixed large sample plot is arranged according to the coordinates of the control sample points.

2. The method for laying out a large forest fixed plot using RTK according to claim 1, characterized in that: The method for filtering out noise points in point cloud data is based on a statistical filtering algorithm. The algorithm formula is as follows: i (x i ,y i ,z i ), calculate the average distance to the surrounding k neighboring points And the average distance of all points is μ d and standard deviation σ d ; like Then the decision point P i The noise points are filtered out, where n is the set threshold coefficient, which is 2-3; The method for removing isolated points from the point cloud data is based on a clustering analysis algorithm. For any point P in the point cloud data i , calculate the number of neighboring points N(P i ,r) within a radius r. If N(P i ,r) < m, then mark the point P i as an isolated point and remove it, where m is the minimum number of cluster points threshold.

3. The method for laying out a large forest fixed plot using RTK according to claim 2, characterized in that: In S1, the UAV is a fixed-wing UAV, and the digital elevation model generation method is based on irregular triangulated network TIN. The generation method is as follows: first, the processed point cloud data is subjected to Delaunay triangulation to construct an irregular triangulated network. Assume that any triangle T in the triangulated network abc , 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) in the triangle, its elevation z p Calculated by linear interpolation: Generate a digital elevation model by traversing the entire triangulated area; The digital ground points are 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 in the area, and the digital height of the digital ground model is added to the digital ground points in the thinned digital ground model.

4. The method for laying out a large forest fixed plot using RTK according to claim 3, characterized in that: Methods for planning a field plot for the three-dimensional digital model of the original terrain and calculating the center of the UAV's outward expansion route through the center of the field plot include: a1. Based on the three-dimensional digital model, the layout area of ​​the fixed large sample plot is divided into a plurality of field sample plots; a2. determining the field sample plot scope of the field sample plot, wherein the field sample plot scope includes the field sample plot area and the field sample plot boundary; a3. Obtain the center coordinates of the field sample plot; a4. Based on the center coordinates of the field sample plot, obtain the center coordinates of the UAV's outward expansion route.

5. The method for laying out a large forest fixed plot using RTK according to claim 4, characterized in that: The method of back-calculating the center of the UAV's outward expansion route from the center of the actual sample plot is to back-calculate the center of the fixed large sample plot from the coordinates of the center of the actual sample plot within the actual sample plot range. Suppose the coordinates of the center of the actual sample plot are (X si ,Y si ), where i = 1, 2, ..., m, m is the number of field plots. First, the average value X of the center coordinates of all field plots in the X direction and the average value Y in the Y direction are calculated as follows: Then, As the initial estimate, based on the coordinates of the expanded sample points (X spj ,Y spj ,H spj ), j = 1, 2, ..., n, n is the number of external expansion sample points, and weighted least squares adjustment optimization is performed; Assume that the distance from the expanded sample point to the initial estimated center is Construct the error equation: where v Xj 、v Yj is the observation residual in the X and Y directions, δX and δY are the correction values ​​of the center coordinates in the X and Y directions, Δd j The correction number for the distance from the expanded sample point to the optimized center coordinate; According to the least squares principle Construct the normal equation to solve δX and δY, where p j The weight of the external sample point is determined according to the accuracy of the point, and the normal equation is constructed to solve δX and δY; Finally, the center coordinates of the large sample site (X c ,Y c ) is calculated as: The center coordinates of the fixed large sample plot are superimposed on the three-dimensional digital model of the original terrain to obtain the center coordinates of the outward expansion route; Assume that the center coordinates of the actual sample plot are (X si ,Y si ), the offset of the center of the fixed large plot relative to the center of the actual plot in the X direction is ΔX, and the offset in the Y direction is ΔY, then the coordinates of the center of the UAV's outward expansion route (X o ,Y o ) is calculated as: X o =X si +ΔX,Y o =Y si +ΔY。 6. The method for laying out a large forest fixed plot using RTK according to claim 5, characterized in that: The method for determining the outward expansion flight altitude of the UAV through the outward expansion route center 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 outward expansion route center according to the height difference, and set the flight altitude to the elevation value from the lowest point in the fixed large sample plot to the outward expansion route center according to the distance; Let the coordinates of the center of the actual plot be , and let the elevation of the highest point in the fixed large plot be H max , the lowest point elevation is H min , the horizontal distance from the lowest point in the plot to the center of the outward expansion route is D, then the outward expansion flight altitude H filght The calculation formula is: Where (ΔH) = H max -H min .

7. The method for laying out a large forest fixed plot using RTK according to claim 6, characterized in that: The method for calculating the coordinates and elevation of the expanded sample points through the expanded sample points includes: b1. Using the center coordinates of the fixed large sample plot to overlay the three-dimensional digital model of the original terrain, aligning the center coordinates of the fixed large sample plot with the outward expansion flight route to obtain the outward expansion flight route; b2. Use a drone to fly at the altitude value according to the outward expansion flight route to perform outward expansion scanning and obtain the coordinates and elevation of the outward expansion sample points; Assume that the center coordinates of the extended flight route are (X o ,Y o ), the offset distance of the expanded sample point relative to the center of the expanded flight route in the X direction is ΔX spj , the offset distance in the Y direction is ΔY spj , the radius of the earth is R, the angle between the line connecting the outward expansion point and the center of the outward expansion flight route and the east direction is θ, then the coordinates of the outward expansion point (X spj ,Y spj ) is calculated as: 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 According to the set expansion distance interval and the expansion point sequence number, let the expansion step length be d, and the sequence number of the nth expansion point in this direction be n, then ΔX spj =n×d×cosα,ΔY spj =n×d×sinα; Assume that the flight altitude is H flight , the center elevation of the extended flight route is H o The vertical height deviation of the expanded sample point relative to the center of the expanded flight route is ΔH spi , then the elevation of the extended sample point H spi The calculation formula is: spi =H o +H flight +ΔH spi .

8. The method for laying out a large forest fixed plot using RTK according to claim 7, characterized in that: The outward expansion flight route is expanded outward from the center of the outward expansion flight route, and the scanning point is the outward expansion sample point. The method for calculating the coordinates and elevation of the outward expansion sample point through the outward expansion sample point is to obtain a 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; 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: In latitude Direction: In the longitude λ direction: 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 It is the center longitude of the outward flight route.

9. The method for laying out a large forest fixed plot using RTK according to claim 8, characterized in that: The coordinates and elevations of the expanded sample points are calculated by inferring the coordinates of the control sample points of the fixed large sample site. The calculation method is as follows: 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: Solve for a, b, c; Assume that the fixed plot is a rectangle with length L and width W, and fit the coordinates of the center point of the plot on the plane (X c ,Y c ,H c ) is used as the reference, and m control sample points are set at equal intervals in the long direction i. In the width direction j, k control sample points are set at equal intervals. Then the coordinates of the control sample point in the jth row and the ith column (X cpij ,Y cpij ,H cpij ) is calculated as follows: X cpij =X c +(i-1)ΔX, where (i=1,2,…,m) Y cpij =Y c +(j-1)ΔY, where (j=1,2,…,k) H cpij =aX cpij +bY cpij +c The pole markers are set at the control sample points by using the GNSS instrument to obtain the coordinates of the pole markers, and the layout direction of the pole markers is set to be consistent with the fixed sample plot.

10. The method for laying out a large forest fixed plot using RTK according to claim 9, characterized in that: The method also includes calculating the accuracy of the coordinates of the pole markers, and calculating the accuracy of the coordinates of the pole markers and the coordinates of the control sample points through the external expansion sample points and the external expansion flight route. The calculation method is: assuming that the theoretical coordinates of the control sample points calculated through the external expansion sample points are (X ct ,Y ct ,H ct ), the actual pole coordinate is (X rb ,Y rb ,H rb ), the precision threshold is ε, then the precision value P is calculated as: When P<ε, the pole layout meets the accuracy requirements; when P>ε, the pole layout does not meet the accuracy requirements.

Images