A terrain surveying and mapping system and a terrain surveying and mapping method based on aerial surveying and mapping technology

Abstract

The present invention provides a terrain surveying and mapping system and a surveying and mapping method based on aerial surveying and mapping technology, which belongs to the technical field of land surveying and mapping. The terrain surveying and mapping method includes: obtaining multiple image data and radar point cloud data of the area to be surveyed, preprocessing the image data of the area to be surveyed, and each image data of the area to be surveyed includes position information and elevation information; constructing a convolutional neural network model for surveying and mapping terrain; obtaining the weather index of the area to be surveyed; obtaining the fitting curve of the flight altitude and terrain complexity, weather index, and camera parameters through experiments; the unmanned aerial vehicle flies according to the determined flight altitude and position according to one-to-one matching, and obtains the image of the new area to be surveyed; the image of the new area to be surveyed is compared, filled, and combined with the image of the area to be surveyed to obtain a surveying image; importing the obtained radar point cloud data; performing aerial triangulation on the surveying image, and establishing a three-dimensional terrain map in combination with the three-dimensional contour of the surface.

Description

Technical Field

[0001] The present application relates to the field of land surveying and mapping technology, and specifically to a terrain surveying and mapping system and method based on aerial surveying and mapping technology. Background Art

[0002] The topographic mapping of a city is the basic technical data necessary for urban planning and development, municipal engineering design, and is also a mapping project with a large capital investment. The original mapping was mainly simulated by optical and mechanical instruments such as theodolites, levels, and flatbeds, combined with manual records, field mapping, and indoor mapping. With the application of global satellite positioning systems and the development of computer technology and network technology, mapping technology is also developing in the direction of digitalization, networking, and convenience. There are more and more mapping applications that integrate airborne LiDAR, camera, and satellite positioning. However, in the current information-based mapping, on the one hand, when stitching images, it is often necessary for the overlap rate of two images to be more than 60%, or to use computer algorithms to extract features from images and then stitch them according to the overlap of coordinates of the same features. On the other hand, when the mapping range is large, factors such as different mapping heights and camera viewing angles affect the actual mapping accuracy, and the amount of data generated is large, which puts a lot of pressure on data transmission and processing. Therefore, how to simplify the mapping workload and ensure that the obtained images maintain a higher accuracy is a question worth considering.

[0003] The patent with the prior art publication number CN114562977A provides a method for urban terrain mapping based on the Internet of Things. It uses actual terrain factors such as roads and paths in the mountains and forests to select convenient reference points to map the urban area.

[0004] The plan is divided into multiple small areas, and then the Internet of Things technology is used to carry out topographic mapping of each small area from the sub-area to the overall area. Image processing is convenient, equipment cost is low, and field mapping work is simple. Moreover, based on drone operation, short-distance and small-scale flight mapping does not require the use of professional aviation equipment such as aircraft, which greatly reduces the cost of mapping and the convenience of operation. Moreover, the mapping map of each small area can adjust the mapping route and the amount of mapping data in real time according to the actual terrain complexity, which is conducive to improving the accuracy of the final mapping.

[0005] Although the above-mentioned existing technical solutions divide the urban plan into multiple small areas, and then use the Internet of Things technology to map the terrain of each small area based on drone operations, the following problems still exist: when mapping each small area, since there is no manual on-site intervention to adjust the flight path of the drone, the accuracy of the detailed representation of the complex terrain in the survey map may be reduced. Summary of the invention

[0006] In order to overcome a series of defects in the prior art, the purpose of this patent is to provide a terrain mapping method based on aerial mapping technology to address the above problems, including:

[0007] Step S1, obtaining a plurality of image data and radar point cloud data of the area to be surveyed and mapped, and preprocessing the image data of the area to be surveyed and mapped, wherein each image data of the area to be surveyed and mapped includes location information and elevation information;

[0008] Step S2, constructing a convolutional neural network model for mapping terrain, optimizing the convolutional neural network model, inputting the preprocessed images of the area to be mapped into the optimized convolutional neural network model in sequence, extracting the features of the preprocessed images of the area to be mapped, predicting the terrain complexity of the images of the area to be mapped based on the extracted features, and obtaining the terrain complexity of the multiple areas to be mapped;

[0009] Step S3, obtaining weather index of the area to be surveyed, including wind speed, visibility, and light intensity; obtaining camera parameters of the aerial surveying camera, including focal length, pixel, aperture, and sensitivity;

[0010] Step S4, obtaining a fitting curve of flight altitude, terrain complexity, weather index, and camera parameters through experiments, substituting the obtained terrain complexity and weather index into the fitting curve, obtaining multiple predetermined flight altitudes of the target area, comparing the predetermined flight altitude with the elevation information of the target area, if the predetermined flight altitude is greater than the altitude of the highest point of the target area, then using the predetermined flight altitude as the determined flight altitude, if the predetermined flight altitude is less than the altitude of the highest point of the target area, then using the altitude of the highest point of the target area + 5m as the determined flight altitude, and matching the determined flight altitude of the target area with the location information;

[0011] Step S5, the UAV flies according to the one-to-one matching determined flight altitude and position to obtain an image of the new area to be surveyed;

[0012] Step S6, comparing, filling, and combining the new image of the area to be surveyed and mapped with the image of the area to be surveyed and mapped obtained in step S1 to obtain a surveying image;

[0013] Step S7, import the acquired radar point cloud data, convert the point_cloud variable type from Numpy to Open3Do3d.geometry.PointCloud type for further processing, select the gridding strategy for surface reconstruction, add a clipping step to remove the unwanted artifacts highlighted in yellow in the figure in order to obtain a clean result, and preliminarily establish the three-dimensional surface contour of the area to be surveyed;

[0014] Step S8, performing aerial triangulation on the surveying and mapping image, and establishing a three-dimensional topographic map in combination with the three-dimensional contour of the surface.

[0015] Furthermore, in step S1, the method for obtaining multiple image data and radar point cloud data of the area to be surveyed is:

[0016] The drone equipped with photography equipment and aperture radar makes a preliminary flight at a preset flight altitude to obtain multiple image data and radar point cloud data of the area to be surveyed;

[0017] Furthermore, in step S1, the method for preprocessing the image of the area to be surveyed is:

[0018] The image of the area to be surveyed is converted into a grayscale image using the weighted average method, the grayscale value is interpolated using the nearest neighbor interpolation algorithm, and the image noise is removed and the image is sharpened using a Gaussian filter to highlight the edge contour of the terrain and obtain a processed image;

[0019] Furthermore, in step S2, the method for constructing a convolutional neural network model for mapping terrain includes the following steps:

[0020] Step S2.1, use drones to collect images of different terrains, and calculate the terrain complexity in the images through field measurements. The terrain complexity is calculated mainly by calculating four indicators: local height difference, local standard deviation, local wrinkle degree, and local total curvature. Among them, the local height difference calculates the difference between the maximum and minimum values ​​of the local window elevation, reflecting the degree of undulation of the local window terrain; the local standard deviation calculates the standard deviation of the local window elevation of the terrain, reflecting the average degree of undulation of the local window terrain; the local wrinkle degree calculates the ratio of the three-dimensional surface area of ​​the local terrain to the projection plane area, reflecting the average degree of wrinkle of the local window terrain; the local total curvature describes the average curvature of the surface eye in all directions, reflecting the average mutation degree of the local window terrain, and its expression is: c=(r2+2s2+t2)1 / 2, where r represents the second-order partial derivative in the x direction, t represents the second-order partial derivative in the y direction, and s represents the second-order partial derivative in the x and y directions: r=fxx=∂2f / ∂x2,t=fyy=∂2f / ∂y2,s=fxy=∂2f / ∂xy; The calculated local height difference, local standard deviation, local wrinkle degree, and local total curvature are normalized, and the formula is: N=(T-minT) / (maxT-minT), T is the original local height difference, local standard deviation, local wrinkle degree, and local total curvature, minT is the minimum local height difference, local standard deviation, local wrinkle degree, and local total curvature value, maxT is the maximum local height difference, local standard deviation, local wrinkle degree, and local total curvature value. After normalization, the local height difference NLRE, local standard deviation NLSD, local wrinkle degree NLRU, and local total curvature NLTC are obtained, and the terrain complexity C=(NLRE+NLSD+NLRU+NLTC) / 4. The collected images are manually annotated with terrain complexity. The images are used as features and the terrain complexity is used as labels to construct a data set, which is divided into a training set and a validation set.

[0021] Step S2.2, build a convolutional neural network architecture, based on the classic convolutional neural network model VGG-16, add convolutional layers, pooling layers, fully connected layers and regression layers, calculate the model's mean absolute error MAE and mean relative error MRE as the loss function, the calculation formula is: MAE=∑mi=1|Mi-Mj|, MRE=∑mi=1(|Mi-Mj| / Mi)*100%, Mi represents the measured value of terrain complexity (%), Mj represents the predicted value of terrain complexity; input the training set into the convolutional neural model for training, and use the back propagation algorithm to update the weights of the convolutional neural network; use the validation set to optimize the hyperparameters of the convolutional neural network, including the learning rate, number of iterations, and optimizer, and obtain the optimal convolutional neural network model with the target task of predicting terrain complexity;

[0022] Furthermore, the method of optimizing the hyperparameters of the convolutional neural network using the validation set includes automatic parameter tuning using the KerasTunner automatic parameter tuning tool and manual parameter tuning by changing the optimizer. Keras Tunner is a library that can automatically search for model training parameters. Its basic idea is to insert a special object (with a specified parameter range) where the parameters need to be adjusted, and then call a search method similar to training.

[0023] Furthermore, in step S2, the method for obtaining the terrain complexity of the plurality of areas to be surveyed is: inputting the plurality of images of the areas to be surveyed into the optimized convolutional neural network model in sequence, the convolutional neural network model automatically extracts the features of the preprocessed images of the areas to be surveyed, and predicts the terrain complexity of the images of the areas to be surveyed based on the extracted features, thereby obtaining the terrain complexity of the plurality of areas to be surveyed;

[0024] Furthermore, in step S4, the method for obtaining the fitting curve of flight altitude, terrain complexity, weather index, and camera parameters through experiments is as follows: the wind speed, visibility, light intensity, and camera in the weather index are taken as variables, and the influence of the weather index on the surveying and mapping accuracy is studied by the control variable method to obtain the fitting curve of the weather index: fweather index=ω11x+ω12y+ω13z, where x is the wind speed, y is the visibility, z is the light intensity, ω11 is the weight of the wind speed, ω12 is the weight of the visibility, and ω13 is the weight of the light intensity; the focal length, pixel, aperture, and sensitivity in the camera parameters are taken as variables, and the influence of the camera parameters on the surveying and mapping accuracy is studied by the control variable method. , obtain the fitting curve of camera parameters fcamera parameters = ω21a + ω22b + ω23c + ω24d, a is focal length, b is pixel, c is aperture, d is sensitivity, ω21 is the weight of focal length, ω22 is the weight of pixel, ω23 is the weight of aperture, ω24 is the weight of sensitivity; take flight altitude, terrain complexity, weather index, and camera parameters as variables, and obtain the fitting curve of flight altitude, terrain complexity, weather index, and camera parameters through the control variable method f = ω1f weather index + ω2f camera parameters + ω3f terrain complexity, ω1 is the weight of weather index, ω2 is the weight of camera parameters, ω3 is the weight of terrain complexity, and f is flight altitude;

[0025] Furthermore, the method for obtaining the fitting curve is as follows: assuming that the data set is {(x1, y1), (x2, y2), (xn, yn)}, xi = (xi1, xi2, xi3, ..., xim), i∈{1, 2, 3, ..., n}, n is the number of data sets, m is the number of variables, and the least squares parameter estimation is used to minimize E(ω, b) = ∑ mi = 1 (yi - ωxi - b) 2, ω = (ω1, ω2, ω3, ..., ωm), and E (ω,b) is differentiated with respect to ω and b respectively, and we get E(ω,b)ω=2(ω∑mi=1xi2-∑mi=1(yi-b)xi), E(ω,b)b=2(mb-∑mi=1(yi-ωxi)), and let the partial derivatives be 0, and we get ω=∑mi=1yi(xi-x) / ∑mi=1xi2-1 / m*∑mi=1xi2)2, b=1 / m*∑mi=1(yi-ωxi), where x=1 / m*∑mi=1xi;

[0026] Furthermore, in step S6, the method of comparing, filling and combining the new image of the area to be surveyed and mapped with the image of the area to be surveyed and mapped obtained in step S1 to obtain a surveying and mapping image includes the following steps:

[0027] Step S6.1, compare the image size of the new area to be surveyed with the image size of the area to be surveyed obtained in step S1, and then scale the image of the new area to be surveyed, the formula is: preset flight altitude / determined flight altitude=image size of the area to be surveyed obtained in step S1 / image size of the new area to be surveyed;

[0028] Step S6.2, matching the feature points of the new image size of the area to be surveyed with the image of the area to be surveyed obtained in step S1, extracting feature points and descriptors from the image using SIFT scale-invariant feature transform, matching the feature points of the two images using brute force matching or nearest neighbor search, removing false matches using a distance threshold-based screening method, filtering duplicate and fuzzy matches using a post-processing method based on the RANSAC algorithm, and obtaining matching results;

[0029] Step S6.3, cutting out the image at the corresponding position of the image of the area to be surveyed and mapped obtained in step S1 according to the matching result, filling the scaled new image of the area to be surveyed and mapped into the position after the image of the area to be surveyed and mapped obtained in step S1 after cutting by image filling, and performing color and space correction;

[0030] Furthermore, aerial triangulation includes: determining relative orientation results according to digital aerial images; performing absolute orientation in combination with IMU / DGPS / GPS data and image control point coordinates; and calculating aerial triangulation results based on the absolute orientation results.

[0031] Furthermore, during aerial triangulation, DEM is collected; DOM is produced based on the DEM collection results; TDOM is produced based on DOM and ground elevation vector data; DEM and TDOM are superimposed on the three-dimensional contour of the surface to construct a three-dimensional topographic map.

[0032] Furthermore, when establishing the three-dimensional contour of the surface, DEM data is also collected.

[0033] The present application also provides a terrain surveying and mapping system based on aerial surveying and mapping technology, including:

[0034] Data acquisition module, used to collect image data and radar point cloud data of the area to be surveyed;

[0035] An image preprocessing module is used to preprocess the collected image data of the area to be surveyed and mapped;

[0036] A deep learning module for predicting terrain complexity based on images of the area to be mapped;

[0037] The information collection module is used to obtain the weather index of the area to be surveyed and mapped, as well as the camera parameters of the aerial surveying and mapping camera. The weather index includes wind speed, visibility, and light intensity, and the camera parameters include focal length, pixels, aperture, and sensitivity.

[0038] A flight altitude adjustment module, used to adjust the flight altitude according to the acquired information and data;

[0039] An image processing module for comparing, filling and combining the images obtained from the first flight and the second flight;

[0040] The 3D modeling module is used to create a 3D topographic map of the area to be surveyed.

[0041] Compared with the prior art, this application has at least the following technical effects or advantages:

[0042] 1) The present invention uses a convolutional neural network based on deep learning to construct a convolutional neural network model with the target task of predicting terrain complexity on the basis of the classic convolutional neural network model VGG-16. The convolutional neural network model inputs the image of the area to be surveyed and mapped, and outputs the terrain complexity of the surveyed area, which eliminates the process of manual field measurement and calculation, shortens the working time of surveying and mapping, and improves the working efficiency of surveying and mapping.

[0043] 2) The present invention takes into account the influence of terrain complexity, weather index and the camera itself on the surveying accuracy, and through a series of analytical calculations, conducts two surveys, automatically adjusts the flight altitude of the drone during the second survey, and compares, fills and combines the images obtained from the two surveys, thereby improving the accuracy of the surveying work.

[0044] 3) The present invention utilizes photographic equipment and aperture radar to simultaneously acquire two data sources, namely, image data of the area to be surveyed and mapping and radar point cloud data, in one flight. Multiple field work and surveying are not required. The two data can be processed simultaneously, which can improve the utilization rate of surveying and mapping data. It has obvious advantages in production efficiency and production cost, and can also avoid consuming a lot of manpower and material resources and reduce measurement costs. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] Figure 1 A module diagram of a terrain surveying and mapping system based on aerial surveying and mapping technology according to the present invention;

[0046] Figure 2 A flowchart of the steps of a terrain surveying and mapping method based on aerial surveying and mapping technology of the present invention;

[0047] Figure 3 Step S2.1 of the terrain surveying and mapping method based on aerial surveying and mapping technology of the present invention is to collect a terrain complexity statistics table calculated by measurement;

[0048] Figure 4 This is a model parameter adjustment optimization effect diagram of step S2.2 of a terrain surveying and mapping method based on aerial surveying and mapping technology of the present invention;

[0049] Figure 5This is a comparison chart of the effect of manually replacing optimizer parameters in a terrain surveying and mapping method based on aerial surveying and mapping technology of the present invention. DETAILED DESCRIPTION

[0050] In order to make the purpose, technical scheme and advantages of the implementation of the present invention clearer, the technical scheme in the embodiment of the present invention will be described in more detail below in conjunction with the drawings in the embodiment of the present invention. In the drawings, the same or similar reference numerals throughout represent the same or similar elements or elements with the same or similar functions. The described embodiments are part of the embodiments of the present invention, not all of the embodiments.

[0051] Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in the field without making any creative work shall fall within the scope of protection of the present invention.

[0052] The embodiments and directional terms described below with reference to the accompanying drawings are exemplary and intended to be used to explain the present invention, but should not be construed as limiting the present invention.

[0053] A terrain surveying and mapping method based on aerial surveying and mapping technology, comprising:

[0054] In this embodiment, a certain surveying area is surveyed, the survey subject is the terrain of the entire area, the target area is surveyed at a scale of 1:1000, and the size of the surveying area is 25.3km. Figure 2 As shown, the following steps are included:

[0055] Step S1, obtaining a plurality of image data and radar point cloud data of the area to be surveyed and mapped, and preprocessing the image data of the area to be surveyed and mapped, wherein each image data of the area to be surveyed and mapped includes location information and elevation information;

[0056] Step S2, constructing a convolutional neural network model for mapping terrain, optimizing the convolutional neural network model, inputting the preprocessed images of the area to be mapped into the optimized convolutional neural network model in sequence, extracting the features of the preprocessed images of the area to be mapped, predicting the terrain complexity of the images of the area to be mapped based on the extracted features, and obtaining the terrain complexity of the multiple areas to be mapped;

[0057] Step S3, obtaining weather index of the area to be surveyed, including wind speed, visibility, and light intensity; obtaining camera parameters of the aerial surveying camera, including focal length, pixel, aperture, and sensitivity;

[0058] Step S4, obtaining a fitting curve of flight altitude, terrain complexity, weather index, and camera parameters through experiments, substituting the obtained terrain complexity and weather index into the fitting curve, obtaining multiple predetermined flight altitudes of the target area, comparing the predetermined flight altitude with the elevation information of the target area, if the predetermined flight altitude is greater than the altitude of the highest point of the target area, then using the predetermined flight altitude as the determined flight altitude, if the predetermined flight altitude is less than the altitude of the highest point of the target area, then using the altitude of the highest point of the target area + 5m as the determined flight altitude, and matching the determined flight altitude of the target area with the location information;

[0059] Step S5, the UAV flies according to the one-to-one matching determined flight altitude and position to obtain an image of the new area to be surveyed;

[0060] Step S6, comparing, filling, and combining the new image of the area to be surveyed and mapped with the image of the area to be surveyed and mapped obtained in step S1 to obtain a surveying image;

[0061] Step S7, import the acquired radar point cloud data, convert the point_cloud variable type from Numpy to Open3Do3d.geometry.PointCloud type for further processing, select the gridding strategy for surface reconstruction, add a clipping step to remove the unwanted artifacts highlighted in yellow in the figure in order to obtain a clean result, and preliminarily establish the three-dimensional surface contour of the area to be surveyed;

[0062] In this example, the meshing strategy using Poisson reconstruction is a bit more technical / mathematical. Its approach is called implicit meshing method, which can be described as trying to "wrap" the data in a smooth cloth. Without going into too much detail, it tries to fit a watertight surface to the original point set by creating a brand new point set representing an isosurface linked with normals. To get the result using Poisson, you just need to adjust the parameters passed to the function.

[0063] Step S8, performing aerial triangulation on the surveying and mapping image, and establishing a three-dimensional topographic map in combination with the three-dimensional contour of the surface.

[0064] In this embodiment, the method of obtaining multiple image data and radar point cloud data of the area to be surveyed and mapped in step S1 includes: a drone equipped with a photographic device and an aperture radar pre-flies a trip at a preset flight altitude to obtain multiple image data and radar point cloud data of the area to be surveyed and mapped;

[0065] In actual surveying, the relationship between flight altitude and accuracy is approximately linear. When the reflectivity of the ground is 10%, the detection distance is 190 meters. When the altitude accuracy is guaranteed, the recommended flight altitude is set to 50-100 meters, and the recommended flight altitude is not more than 150 meters. In this embodiment, the flight altitude is pre-set to 100 meters, the heading overlap is 75%, the side phase overlap is 55%, and the flight speed is 15m / s.

[0066] In this embodiment, the method for preprocessing the image of the area to be surveyed and mapped in step S1 includes: using a weighted average method to convert the image of the area to be surveyed and mapped into a grayscale image, using a nearest neighbor interpolation algorithm to interpolate grayscale values, using a low-pass filter such as a Gaussian filter to remove image noise and perform image sharpening, highlighting the edge contour of the terrain, and obtaining a processed image;

[0067] In this embodiment, a total of 569 pieces of image data are acquired in the first flight, and grayscale conversion, image denoising and image sharpening are performed on the acquired image data.

[0068] In this embodiment, the method of constructing a convolutional neural network model for surveying and mapping terrain in step S2 includes:

[0069] Step S2.1, use drones to collect images of different terrains, and calculate the terrain complexity in the images through field measurements. The terrain complexity is calculated mainly by calculating four indicators: local height difference, local standard deviation, local wrinkle degree, and local total curvature. Among them, the local height difference calculates the difference between the maximum and minimum values ​​of the local window elevation, reflecting the degree of undulation of the local window terrain; the local standard deviation calculates the standard deviation of the local window elevation of the terrain, reflecting the average degree of undulation of the local window terrain; the local wrinkle degree calculates the ratio of the three-dimensional surface area of ​​the local terrain to the projection plane area, reflecting the average degree of wrinkle of the local window terrain; the local total curvature describes the average curvature of the surface eye in all directions, reflecting the average mutation degree of the local window terrain, and its expression is: c=(r2+2s2+t2)1 / 2, where r represents the second-order partial derivative in the x direction, t represents the second-order partial derivative in the y direction, and s represents the second-order partial derivative in the x and y directions: r=fxx=∂2f / ∂x2,t=fyy=∂2f / ∂y2,s=fxy=∂2f / ∂xy; The calculated local height difference, local standard deviation, local wrinkle degree, and local total curvature are normalized, and the formula is: N=(T-minT) / (maxT-minT), T is the original local height difference, local standard deviation, local wrinkle degree, and local total curvature, minT is the minimum local height difference, local standard deviation, local wrinkle degree, and local total curvature value, maxT is the maximum local height difference, local standard deviation, local wrinkle degree, and local total curvature value. After normalization, the local height difference NLRE, local standard deviation NLSD, local wrinkle degree NLRU, and local total curvature NLTC are obtained, and the terrain complexity C=(NLRE+NLSD+NLRU+NLTC) / 4. The collected images are manually annotated with terrain complexity. The images are used as features and the terrain complexity is used as labels to construct a data set, which is divided into a training set and a validation set.

[0070] In this embodiment, the terrain complexity statistics table calculated by measurement is collected as follows: Figure 3 shown.

[0071] In this embodiment, the Hobson method is used to calculate the local wrinkle degree. The center points of 9 pixels form 8 triangles. The sum of the areas of the triangles represents the surface area As of the area formed by the center points of the 9 pixels. The horizontal projection area is Ap (the area of ​​4 grids), so the wrinkle degree is: TLRU=As / Ap=As / 4g2.

[0072] In this embodiment, the key to calculating the full area rate is to determine the second-order derivative and partial derivative of the central pixel. In digital terrain analysis, the derivative on the grid DEM can be approximated by numerical difference, that is:

[0073] fx=(H3+H6+H9-H1-H4-H7) / 6g

[0074] fy=(H7+H8+H9-H1-H2-H3) / 6g

[0075] fxx=(H1+H3+H4+H6+H7+H9-2(H2+H5+H8)) / 3g2

[0076] fyy=(H1+H2+H3+H7+H8+H9-2(H4+H5+H8)) / 3g2

[0077] fxy=fyx=(H9+H1-H7-H3) / 4g2

[0078] Where H1, H2, and H9 represent the pixel heights of the local window.

[0079] In this embodiment, a total of 564 images of different terrain complexity are collected and constructed into a data set after manual annotation, 80% of which are divided into a training set and 20% are divided into a validation set.

[0080] Step S2.2, build a convolutional neural network architecture, based on the classic convolutional neural network model VGG-16, add convolutional layers, pooling layers, fully connected layers and regression layers, calculate the model's mean absolute error MAE and mean relative error MRE as the loss function, the calculation formula is: MAE=∑mi=1|Mi-Mj|, MRE=∑mi=1(|Mi-Mj| / Mi)*100%, Mi represents the measured value of terrain complexity (%), Mj represents the predicted value of terrain complexity; input the training set into the convolutional neural model for training, and use the back propagation algorithm to update the weights of the convolutional neural network; use the validation set to optimize the hyperparameters of the convolutional neural network, including the learning rate, number of iterations, and optimizer, and obtain the optimal convolutional neural network model with the target task of predicting terrain complexity. The model parameter optimization effect diagram is shown in the figure. Figure 4 shown.

[0081] In this embodiment, the terrain complexity R2 predicted by the model is 0.928, the mean absolute error MAE is 6.1%, the mean relative error is 8.9%, the optimized learning rate is 0.001, and the number of batch data is 64.

[0082] In this embodiment, the method of optimizing the hyperparameters of the convolutional neural network using the validation set includes automatic parameter adjustment using the KerasTunner automatic parameter adjustment tool and manual parameter adjustment by the optimizer. Keras Tunner is a library that can automatically search for model training parameters. Its basic idea is to insert a special object (the parameter range can be specified) where the parameters need to be adjusted, and then call a search method similar to training.

[0083] In this embodiment, during the parameter optimization of the model, this embodiment determines whether it is optimal by studying the function loss of the training set and the validation set. In the early stage, there is a certain fluctuation in the model loss. As the number of training epochs increases, the loss of the training set and the validation set gradually decreases, and the average error of both is below 0.1. It can be seen that the best number of training times is about 50.

[0084] In this example, the optimizer is manually adjusted and four optimizers are tested, namely SGD, DOMentum, RMSprop, and Adam. The results are as follows: Figure 5 As shown, it can be seen that the Adam optimizer has the fastest fitting and the best effect.

[0085] In this embodiment, in step S2, the method for obtaining the terrain complexity of multiple areas to be surveyed is: inputting multiple images of the areas to be surveyed into the optimized convolutional neural network model in sequence, the convolutional neural network model automatically extracts the features of the pre-processed images of the areas to be surveyed, and predicts the terrain complexity of the images of the areas to be surveyed based on the extracted features, thereby obtaining the terrain complexity of the multiple areas to be surveyed;

[0086] In this embodiment, in step S4, the method for obtaining the fitting curve of the flight altitude and the terrain complexity, the weather index, and the camera parameters through experiments is as follows: the wind speed, visibility, light intensity, and the camera in the weather index are taken as variables, and the influence of the weather index on the mapping accuracy is studied by the control variable method to obtain the fitting curve of the weather index: fweather index=ω11x+ω12y+ω13z, where x is the wind speed, y is the visibility, z is the light intensity, ω11 is the weight of the wind speed, ω12 is the weight of the visibility, and ω13 is the weight of the light intensity; the focal length, pixel, aperture, and sensitivity in the camera parameters are taken as variables, and the influence of the camera parameters on the mapping accuracy is studied by the control variable method. The fitting curve of camera parameters is obtained: fcamera parameters = ω21a + ω22b + ω23c + ω24d, where a is focal length, b is pixel, c is aperture, d is sensitivity, ω21 is the weight of focal length, ω22 is the weight of pixel, ω23 is the weight of aperture, and ω24 is the weight of sensitivity; the flight altitude, terrain complexity, weather index, and camera parameters are taken as variables, and the fitting curve of flight altitude, terrain complexity, weather index, and camera parameters is obtained by the control variable method: f = ω1fweather index + ω2fcamera parameters + ω3fterrain complexity, where ω1 is the weight of weather index, ω2 is the weight of camera parameters, ω3 is the weight of terrain complexity, and f is the flight altitude;

[0087] In this embodiment, the method for obtaining the fitting curve is: assuming that the data set is {(x1, y1), (x2, y2), (xn, yn)}, xi=(xi1, xi2, xi3, ..., xim), i∈{1, 2, 3, ..., n}, n is the number of data sets, m is the number of variables, and the least squares parameter estimation is used to minimize E(ω, b)=∑mi=1(yi-ωxi-b)2, ω=(ω1, ω2, ω3, ..., ωm), and E(ω,b) is differentiated with respect to ω and b respectively, and we get E(ω,b)ω=2(ω∑mi=1xi2-∑mi=1(yi-b)xi), E(ω,b)b=2(mb-∑mi=1(yi-ωxi)), and let the partial derivatives be 0, and we get ω=∑mi=1yi(xi-x) / ∑mi=1xi2-1 / m*∑mi=1xi2)2, b=1 / m*∑mi=1(yi-ωxi), where x=1 / m*∑mi=1xi;

[0088] In this embodiment, the weight of the wind speed of the weather index fitting curve is obtained to be ω11=0.223, the weight of the visibility is obtained to be ω12=0.324, and the weight of the light intensity is obtained to be ω13=0.453.

[0089] In this embodiment, the camera parameter fitting curve is obtained as follows: the focal length weight ω21=0.13, the pixel weight ω22=0.314, the aperture weight ω23=0.428, and the sensitivity weight ω24=0.128.

[0090] In this embodiment, fitting curves of flight altitude, terrain complexity, weather index, and camera parameters are obtained, with the weight of weather index ω1=0.477, the weight of camera parameters ω2=0.365, and the weight of terrain complexity ω3=0.158.

[0091] In this embodiment, in step S6, the method of comparing, filling, and combining the new image of the area to be surveyed and mapped with the image of the area to be surveyed and mapped obtained in step S1 to obtain a surveying and mapping image includes the following steps:

[0092] Step S6.1, compare the image size of the new area to be surveyed with the image size of the area to be surveyed obtained in step S1, and then scale the image of the new area to be surveyed, the formula is: preset flight altitude / determined flight altitude=image size of the area to be surveyed obtained in step S1 / image size of the new area to be surveyed;

[0093] Step S6.2, matching the feature points of the new image size of the area to be surveyed with the image of the area to be surveyed obtained in step S1, extracting feature points and descriptors from the image using SIFT scale-invariant feature transform, matching the feature points of the two images using brute force matching or nearest neighbor search, removing false matches using a distance threshold-based screening method, filtering duplicate and fuzzy matches using a post-processing method based on the RANSAC algorithm, and obtaining matching results;

[0094] Step S6.3, cutting out the image at the corresponding position of the image of the area to be surveyed and mapped obtained in step S1 according to the matching result, filling the scaled new image of the area to be surveyed and mapped into the position after the image of the area to be surveyed and mapped obtained in step S1 after cutting by image filling, and performing color and space correction;

[0095] In this embodiment, aerial triangulation includes: determining a relative orientation result according to digital aerial images; then performing absolute orientation in combination with IMU / DGPS / GPS data and image control point coordinates; and calculating an aerial triangulation result based on the absolute orientation result.

[0096] In this embodiment, it also includes collecting DEM during aerial triangulation; producing DOM based on the DEM collection results; producing TDOM according to DOM and ground elevation vector data; and superimposing DEM and TDOM on the three-dimensional contour of the surface to construct a three-dimensional topographic map.

[0097] In this embodiment, when establishing the three-dimensional contour of the ground surface, DEM data is also collected.

[0098] On the other hand, Figure 1 As shown, this embodiment also provides a terrain surveying and mapping system based on aerial surveying and mapping technology, including:

[0099] The data acquisition module 1 is used to collect image data and radar point cloud data of the area to be surveyed; the image preprocessing module 2 is used to preprocess the collected image data of the area to be surveyed; the deep learning module 3 is used to predict the terrain complexity based on the image of the area to be surveyed; the information acquisition module 4 is used to obtain the weather index of the area to be surveyed, including the wind speed, visibility, light intensity at this time, and obtain the camera parameters of the aerial surveying camera, including focal length, pixel, aperture, sensitivity; the flight altitude adjustment module 5 is used to adjust the flight altitude according to the acquired information and data; the image processing module 6 is used to compare, fill and combine the images obtained from the first flight and the second flight; the three-dimensional modeling module 7 is used to establish a three-dimensional terrain map of the area to be surveyed.

Claims

1. A terrain surveying and mapping method based on aerial surveying and mapping technology, characterized in that: include: Step S1, obtaining a plurality of image data and radar point cloud data of the area to be surveyed and mapped, and preprocessing the image data of the area to be surveyed and mapped, wherein each image data of the area to be surveyed and mapped includes location information and elevation information; Step S2, constructing a convolutional neural network model for mapping terrain, optimizing the convolutional neural network model, inputting the preprocessed images of the area to be mapped into the optimized convolutional neural network model in sequence, extracting the features of the preprocessed images of the area to be mapped, predicting the terrain complexity of the images of the area to be mapped based on the extracted features, and obtaining the terrain complexity of the multiple areas to be mapped; Step S3, obtaining weather index of the area to be surveyed, including wind speed, visibility, and light intensity; obtaining camera parameters of the aerial surveying camera, including focal length, pixel, aperture, and sensitivity; Step S4, obtaining a fitting curve of flight altitude, terrain complexity, weather index, and camera parameters through experiments, substituting the obtained terrain complexity and weather index into the fitting curve, obtaining multiple predetermined flight altitudes of the target area, comparing the predetermined flight altitude with the elevation information of the target area, if the predetermined flight altitude is greater than the altitude of the highest point of the target area, then using the predetermined flight altitude as the determined flight altitude, if the predetermined flight altitude is less than the altitude of the highest point of the target area, then using the altitude of the highest point of the target area + 5m as the determined flight altitude, and matching the determined flight altitude of the target area with the location information; Step S5, the UAV flies according to the one-to-one matching determined flight altitude and position to obtain an image of the new area to be surveyed; Step S6, comparing, filling, and combining the new image of the area to be surveyed and mapped with the image of the area to be surveyed and mapped obtained in step S1 to obtain a surveying image; Step S7, import the acquired radar point cloud data, convert the point_cloud variable type from Numpy to Open3Do3d.geometry.PointCloud type for further processing, select the gridding strategy for surface reconstruction, add a clipping step to remove the unwanted artifacts highlighted in yellow in the figure in order to obtain a clean result, and preliminarily establish the three-dimensional surface contour of the area to be surveyed; Step S8, performing aerial triangulation on the surveying and mapping image, and establishing a three-dimensional topographic map in combination with the three-dimensional contour of the surface.

2. A terrain surveying and mapping method based on aerial surveying and mapping technology according to claim 1, characterized in that: The method for obtaining multiple image data and radar point cloud data of the area to be surveyed is: The drone equipped with photographic equipment and aperture radar will make a preliminary flight at a preset altitude to obtain multiple image data and radar point cloud data of the area to be surveyed.

3. The terrain surveying and mapping method based on aerial surveying and mapping technology according to claim 1 is characterized in that: The method for preprocessing the image of the surveying area is: The weighted average method is used to convert the image of the area to be surveyed into a grayscale image, the nearest neighbor interpolation algorithm is used to interpolate the grayscale values, and the Gaussian filter is used to remove image noise and sharpen the image to highlight the edge contour of the terrain and obtain a processed image.

4. The terrain surveying and mapping method based on aerial surveying and mapping technology according to claim 1 is characterized in that: In step S2, the method for constructing a convolutional neural network model for mapping terrain includes the following steps: Step S2.1, use drones to collect images of different terrains, and calculate the terrain complexity in the images through field measurements. The terrain complexity is calculated mainly by calculating four indicators: local height difference, local standard deviation, local wrinkle degree, and local total curvature. Among them, the local height difference calculates the difference between the maximum and minimum values ​​of the local window elevation, reflecting the degree of undulation of the local window terrain; the local standard deviation calculates the standard deviation of the local window elevation of the terrain, reflecting the average degree of undulation of the local window terrain; the local wrinkle degree calculates the ratio of the three-dimensional surface area of ​​the local terrain to the projection plane area, reflecting the average degree of wrinkle of the local window terrain; the local total curvature describes the average curvature of the surface eye in all directions, reflecting the average mutation degree of the local window terrain, and its expression is: c = (r 2 +2s 2 +t 2 ) 1 / 2 , where r represents the second-order partial derivative in the x direction, t represents the second-order partial derivative in the y direction, and s represents the second-order partial derivative in the x and y directions: r=f xx =∂ 2 f / ∂x 2 ,t=f yy =∂ 2 f / ∂y 2 ,s=f xy =∂ 2 f / ∂xy; the calculated local height difference, local standard deviation, local wrinkle degree, and local total curvature are normalized, and the formula is: N=(T-minT) / (maxT-minT), T is the original local height difference, local standard deviation, local wrinkle degree, and local total curvature, minT is the minimum local height difference, local standard deviation, local wrinkle degree, and local total curvature value, maxT is the maximum local height difference, local standard deviation, local wrinkle degree, and local total curvature value, and the local height difference N is obtained after normalization. LRE , local standard deviation N LSD , local wrinkle degree N LRU , local total curvature N LTC， Terrain complexity C = (N LRE +N LSD +N LRU +N LTC ) / 4, manually annotate the terrain complexity of the collected images, use the images as features and the terrain complexity as labels to construct a dataset, and divide it into a training set and a validation set; Step S2.2, build a convolutional neural network architecture, based on the classic convolutional neural network model VGG-16, add convolutional layers, pooling layers, fully connected layers, and regression layers, and calculate the model's mean absolute error MAE and mean relative error MRE as the loss function. The calculation formula is: MAE=∑ m i=1 |M i -M j |,MRE=∑ m i=1 (|M i -M j | / M i )*100%,M i represents the measured value of terrain complexity (%), M j Represents the predicted value of terrain complexity; the training set is input into the convolutional neural model for training, and the weights of the convolutional neural network are updated using the back-propagation algorithm; the validation set is used to optimize the hyperparameters of the convolutional neural network, including the learning rate, number of iterations, and optimizer, to obtain the optimal convolutional neural network model with the objective task of predicting terrain complexity.

5. The terrain surveying and mapping method based on aerial surveying and mapping technology according to claim 1, characterized in that: In step S2, the method for obtaining the terrain complexity of multiple areas to be surveyed is: multiple images of the areas to be surveyed are sequentially input into the optimized convolutional neural network model, the convolutional neural network model automatically extracts the features of the preprocessed images of the areas to be surveyed, and the terrain complexity of the images of the areas to be surveyed is predicted based on the extracted features, so as to obtain the terrain complexity of multiple areas to be surveyed.

6. The terrain surveying and mapping method based on aerial surveying and mapping technology according to claim 1, characterized in that: In step S4, the method of obtaining the fitting curve of flight altitude, terrain complexity, weather index, and camera parameters through experiments is as follows: taking the wind speed, visibility, light intensity, and camera in the weather index as variables, studying the influence of the weather index on the mapping accuracy through the control variable method, and obtaining the fitting curve f of the weather index 天气指数 =ω 11 x+ω 12 y+ω 13 z, x is the wind speed, y is the visibility, z is the light intensity, ω 11 is the weight of wind speed, ω 12 is the visibility weight, ω 13 is the weight of light intensity; the focal length, pixel, aperture, and sensitivity in the camera parameters are taken as variables, and the influence of camera parameters on surveying and mapping accuracy is studied by the control variable method, and the fitting curve f of the camera parameters is obtained. 相机参数 =ω 21 a+ω 22 b+ω 23 c+ω 24 d, a is focal length, b is pixel, c is aperture, d is sensitivity, ω 21 is the focal length weight, ω 22 is the pixel weight, ω 23 is the aperture weight, ω 24 is the weight of sensitivity; taking flight altitude, terrain complexity, weather index, and camera parameters as variables, the fitting curve f=ω1f of flight altitude, terrain complexity, weather index, and camera parameters is obtained by the control variable method 天气指数 +ω2f 相机参数 +ω3f 地形复杂度， ω1 is the weight of the weather index, ω2 is the weight of the camera parameters, ω3 is the weight of the terrain complexity, and f is the flight altitude.

7. The terrain surveying and mapping method based on aerial surveying and mapping technology according to claim 1, characterized in that: The method to obtain the fitting curve is: Assume that the data set is {(x1,y 1） ,(x2,y2),(x n ,y n )}, x i =(x i 1 ,x i 2 ,x i 3 ,...,x i m ), i∈{1,2,3,...,n}, n is the number of data sets, m is the number of variables, and the least squares parameter estimation is used so that E(ω,b)=∑ m i=1 (y i -ωx i -b) 2 Minimize, ω = (ω 1 ,ω 2 ,ω 3 ,...,ω m ), taking the derivative of E(ω,b) with respect to ω and b respectively, we get E(ω,b)ω=2(ω∑ m i=1 x i 2 -∑ m i=1 (y i -b)x i ), E(ω,b) b =2(mb-∑ m i=1 (y i -ωx i )), let all partial derivatives be 0, and we get ω=∑ m i=1 y i (x i -x) / ∑ m i=1 x i 2 -1 / m*∑ m i=1 x i 2 ) 2 ,b=1 / m*∑ m i=1 (y i -ωx i ), where x=1 / m*∑ m i=1 x i .

8. The terrain surveying and mapping method based on aerial surveying and mapping technology according to claim 1, characterized in that: In step S6, the method of comparing, filling and combining the new image of the area to be surveyed and mapped with the image of the area to be surveyed and mapped obtained in step S1 to obtain a surveying and mapping image includes the following steps: Step S6.1, compare the image size of the new area to be surveyed with the image size of the area to be surveyed obtained in step S1, and then scale the image of the new area to be surveyed, the formula is: preset flight altitude / determined flight altitude=image size of the area to be surveyed obtained in step S1 / image size of the new area to be surveyed; Step S6.2, matching the feature points of the new image size of the area to be surveyed with the image of the area to be surveyed obtained in step S1, extracting feature points and descriptors from the image using SIFT scale-invariant feature transform, matching the feature points of the two images using brute force matching or nearest neighbor search, removing false matches using a distance threshold-based screening method, filtering duplicate and fuzzy matches using a post-processing method based on the RANSAC algorithm, and obtaining matching results; Step S6.3, according to the matching result, the image of the area to be surveyed and mapped corresponding to the image obtained in step S1 is cut out, and the scaled new image of the area to be surveyed and mapped is filled into the position after the image of the area to be surveyed and mapped obtained in step S1 is cut out by image filling, and color and space correction are performed.

9. The terrain surveying and mapping method based on aerial surveying and mapping technology according to claim 1, characterized in that: Aerial triangulation includes: relative orientation results determined according to digital aerial images; absolute orientation based on IMU / DGPS / GPS data and image control point coordinates; and calculation of aerial triangulation results based on absolute orientation results. During aerial triangulation, DEM is collected; DOM is produced based on the DEM collection results; TDOM is produced based on DOM and ground elevation vector data; DEM and TDOM are superimposed on the three-dimensional contour of the surface to construct a three-dimensional topographic map.

10. A terrain surveying and mapping system based on aerial surveying and mapping technology, characterized in that: include: Data acquisition module, used to collect image data and radar point cloud data of the area to be surveyed; An image preprocessing module is used to preprocess the collected image data of the area to be surveyed and mapped; A deep learning module for predicting terrain complexity based on images of the area to be mapped; The information collection module is used to obtain the weather index of the area to be surveyed and mapped, as well as the camera parameters of the aerial surveying and mapping camera. The weather index includes wind speed, visibility, and light intensity, and the camera parameters include focal length, pixels, aperture, and sensitivity. A flight altitude adjustment module, used to adjust the flight altitude according to the acquired information and data; An image processing module for comparing, filling and combining the images obtained from the first flight and the second flight; The 3D modeling module is used to create a 3D topographic map of the area to be surveyed.

Images