Geological body modeling method and system based on unmanned aerial vehicle laser radar

Abstract

The invention relates to the technical field of photogrammetry, in particular to a geologic body modeling method and system based on an unmanned aerial vehicle laser radar, and the method comprises the following steps: obtaining laser ranging and pulse counting, calling a field programmable gate array to carry out time marking, retrieving adjacent frame data, and generating an instantaneous attitude and coordinates through linear interpolation; a three-dimensional laser point cloud is generated based on installation error matrix conversion, a local gradient is calculated by mapping grids to screen ground points, and a geologic body three-dimensional digital model is constructed by utilizing Delaune triangulation. According to the method, the field programmable gate array is utilized to mark data time, the interpolation algorithm is combined to solve the instantaneous attitude, the three-dimensional point cloud is generated based on the error matrix, the grids are mapped, the gradient is calculated to screen the ground points, the geologic model is constructed through Delaune triangulation, and high-density ground surface data can be obtained without contact measurement. And the modeling reduction degree of the complex terrain is improved.

Description

Technical Field

[0001] This invention relates to the field of photogrammetry technology, and in particular to a geological body modeling method and system based on UAV lidar. Background Technology

[0002] Photogrammetry is a field that utilizes cameras and sensors mounted on aircraft or ground platforms to acquire image data of a target, and then determines the target's three-dimensional shape, location, and size by analyzing the geometric relationships within the images. Traditional geological modeling methods involve surveyors carrying total stations, levels, or handheld GPS devices into the survey area. They manually select points and measure the three-dimensional coordinates of surface features, transmitting the recorded coordinate data to a computer terminal. Drawing software then uses grid interpolation algorithms to generate a digital elevation model or contour map of the terrain.

[0003] Traditional geological body modeling requires surveyors to carry heavy surveying equipment deep into the surveyed area, relying on subjective judgment to select surface feature points one by one for contact measurement. Due to limitations in manpower and terrain accessibility, data collection points are sparse and unevenly distributed. When facing complex and rugged terrain, it is difficult to achieve comprehensive coverage, resulting in the loss of a large amount of terrain detail information. Relying on a limited number of discrete coordinate points to generate a model through grid interpolation makes it difficult to accurately reproduce the surface undulation. Not only is the field operation cycle long and labor-intensive, but the manual reading and data transcription process is also prone to human error, which greatly reduces the accuracy of model construction. It cannot meet the actual engineering needs of high-precision full-element 3D digital reconstruction of large-scale complex terrain. Summary of the Invention

[0004] The purpose of this invention is to address the shortcomings of existing technologies by proposing a geological body modeling method and system based on UAV lidar.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a geological body modeling method based on UAV lidar, comprising the following steps: S1: Acquire laser ranging value, scanning angle, angular velocity, acceleration, position coordinates and pulse count value; read the pulse count value generated by field programmable gate array; perform time-stamping processing on the laser ranging value, the scanning angle, the angular velocity, the acceleration and the position coordinates to generate a time-stamped original measurement dataset; S2: Retrieve the angular velocity, acceleration, and position coordinate data frames that are adjacent to the time of the laser ranging value from the original time-stamped measurement dataset, and use the data frames to perform a linear interpolation algorithm to generate instantaneous attitude and instantaneous three-dimensional spatial coordinates; S3: Calculate the three-dimensional coordinates of the sensor coordinate system based on the laser ranging value and the scanning angle; and generate a three-dimensional laser point cloud dataset by transforming the three-dimensional coordinates of the sensor coordinate system based on the installation error matrix, the instantaneous attitude, and the instantaneous three-dimensional spatial coordinates. S4: Map the three-dimensional laser point cloud dataset to a two-dimensional grid to extract potential ground seed points, calculate the ratio of the elevation difference and horizontal distance between the potential ground seed points and their neighbors to generate local slope values, filter the local slope values ​​based on slope thresholds to determine ground points, and input the ground points into the Delaunay triangulation algorithm to generate a three-dimensional digital model of the geological body.

[0006] As a further aspect of the present invention, the original measurement dataset with timestamps includes the laser ranging value marked with the pulse count value at the same moment, the scanning angle, the angular velocity, the acceleration, and the position coordinates; The instantaneous attitude is a quaternion of the body attitude calculated based on linear interpolation; The instantaneous three-dimensional spatial coordinates include interpolated longitude, latitude, and elevation values; The three-dimensional laser point cloud dataset consists of point cloud data in the world geodetic coordinate system after installation error compensation and coordinate system transformation. The three-dimensional digital model of the geological body includes an irregular triangular mesh generated based on the ground points, triangular patch normal vectors, geological body slope parameters, and slope aspect parameters.

[0007] As a further aspect of the present invention, the process of generating the instantaneous attitude and the instantaneous three-dimensional spatial coordinates specifically includes: Read the laser emission time timestamp corresponding to each laser ranging value in the timestamped raw measurement dataset; Using the laser emission time timestamp as the retrieval reference, a bidirectional search is performed in the time series of the original measurement dataset with timestamps to locate the preceding inertial navigation recording frame at the time before the laser emission time timestamp and the following inertial navigation recording frame at the time after the laser emission time timestamp. Extract the position coordinates, angular velocity and acceleration contained in the preceding inertial navigation recording frame and the following inertial navigation recording frame respectively, and calculate the time difference ratio of the laser emission time timestamp relative to the preceding inertial navigation recording frame; Based on the time difference ratio, a linear position weighting calculation is performed on the position coordinates of the preceding inertial navigation recording frame and the following inertial navigation recording frame to generate the instantaneous three-dimensional spatial coordinates; The body attitude angle is calculated by fusing the angular velocity and the acceleration, and spherical linear interpolation is performed on the attitude quaternions of adjacent moments using the time difference ratio to generate the instantaneous attitude corresponding to the timestamp of the laser emission moment.

[0008] As a further aspect of the present invention, the calculation process of the three-dimensional coordinates of the sensor coordinate system specifically includes: Obtain the laser ranging value and the scanning angle recorded in the original measurement dataset with timestamps, and read the internal optical parameters of the lidar device; The laser ranging value is optically corrected based on the internal optical parameters. The corrected laser ranging value and the scanning angle are decomposed into component data of the laser beam in the rectangular coordinate system defined by the sensor itself using trigonometric function relationships. Based on the horizontal, vertical, and longitudinal components obtained from the decomposition, coordinate vectors are constructed to generate the three-dimensional coordinates of the sensor coordinate system. The calculation of the three-dimensional coordinates of the sensor coordinate system follows the formula below: ; in, This represents the generated three-dimensional coordinate vector of the sensor coordinate system. Represents the horizontal axis coordinate value in the sensor coordinate system. This represents the vertical axis coordinate value in the sensor coordinate system. Represents the vertical axis coordinate value in the sensor coordinate system. This represents the laser ranging value. Represents the scanning angle, This represents the fixed pitch deviation angle of the laser emitter relative to the sensor reference plane.

[0009] As a further aspect of the present invention, the generation process of the three-dimensional laser point cloud dataset specifically includes: Obtain the pre-calibrated installation error matrix, which includes the rotation and translation parameters of the lidar relative to the center of the UAV body; The instantaneous attitude is read and converted into a rotation matrix of the body coordinate system relative to the geodetic coordinate system. The instantaneous three-dimensional spatial coordinates are read as the translation vector of the body center in the geodetic coordinate system. Multiply the three-dimensional coordinates of the sensor coordinate system by the installation error matrix to transform them into intermediate transition coordinates in the body coordinate system. Multiply the intermediate transition coordinates by the rotation matrix and add the translation vector to complete the mapping transformation from the sensor's local space to the world's ground space; The calculation results corresponding to all the laser ranging values ​​are traversed, and all the transformed coordinate points are stacked and stored in chronological order to generate the three-dimensional laser point cloud dataset.

[0010] As a further aspect of the present invention, the extraction process of the potential ground seed points specifically includes: Read the extreme values ​​of the planar coordinates of all point cloud data in the three-dimensional laser point cloud dataset to determine the horizontal rectangular boundary range of the point cloud coverage. Based on preset grid resolution parameters, the horizontal rectangular boundary range is divided into non-overlapping regular two-dimensional grid matrices; Traverse each laser point in the three-dimensional laser point cloud dataset, and map the laser point index to the corresponding two-dimensional grid cell according to the planar coordinates of the laser point; For each two-dimensional grid cell containing the laser point, retrieve the elevation values ​​of all laser points within that two-dimensional grid cell, and filter out the laser point with the minimum elevation value through comparison calculation; The laser point with the minimum elevation value is marked as the representative low point of the two-dimensional grid cell, and the representative low points of all the two-dimensional grid cells are aggregated to generate the potential ground seed point.

[0011] As a further aspect of the present invention, the process of generating the local slope value specifically includes: Construct a KD tree spatial index structure for the potential ground seed points, and for each potential ground seed point, use the KD tree spatial index structure to search for its nearest neighbor set within a set radius. Traverse each neighboring point in the nearest neighbor set and extract the elevation and planar coordinate components of the potential ground seed point and the neighboring point respectively. Calculate the absolute value of the elevation difference between the potential ground seed point and its neighboring point, and calculate the Euclidean distance between the potential ground seed point and its neighboring point projected onto the horizontal plane. The slope tangent is calculated based on the ratio of the absolute value of the elevation difference to the Euclidean distance. The maximum value among the calculation results of all the neighboring points is selected to generate the local slope value. The calculation of the local slope value follows the following formula: ; in, This represents the generated local slope value. The set of indices representing the nearest neighbor set, This represents the elevation value of the potential ground seed point. Representing the The elevation values ​​of the neighboring points. The lateral plane coordinates of the potential ground seed points are represented. The longitudinal plane coordinates of the potential ground seed points are represented. Representing the The lateral plane coordinates of the neighboring points. Representing the The vertical plane coordinate values ​​of the neighboring points.

[0012] As a further aspect of the present invention, the process of determining the ground points specifically includes: Obtain the preset slope threshold, which is set based on the average roughness of the geological type of the survey area; Read the local slope value corresponding to each potential ground seed point, and compare the local slope value with the slope threshold value one by one; If the local slope value is less than or equal to the slope threshold, it is determined that the local terrain of the potential ground seed point is flat, and the attribute identifier of the potential ground seed point is modified to a confirmed ground attribute. If the local slope value is greater than the slope threshold, the potential ground seed point is determined to be a vegetation or non-ground obstacle interference point, and the potential ground seed point is removed from the candidate set. Extract all point data marked with the confirmed ground attribute, perform outlier denoising processing on the point data of the confirmed ground attribute, and generate the ground point.

[0013] As a further aspect of the present invention, the process of generating the three-dimensional digital model of the geological body specifically includes: Project the ground points onto a two-dimensional plane to construct the convex hull boundary of the point set; The Deloitte triangulation algorithm is called to connect the ground points as vertices to generate a non-overlapping triangular topology network that covers the entire convex hull boundary, thus generating an irregular triangular network. Traverse each triangular facet in the irregular triangular network, extract the three-dimensional coordinates of the three vertices, and calculate the normal vector of the triangular facet through vector cross product operation; The slope parameter is calculated based on the vertical component of the normal vector, and the aspect parameter is calculated based on the horizontal projection direction of the normal vector. The irregular triangular mesh, the normal vector of the triangular facet, the slope parameters of the geological body, and the aspect parameters are encapsulated to generate a three-dimensional digital model of the geological body.

[0014] A geological body modeling system based on UAV lidar, the system being used to implement the aforementioned geological body modeling method based on UAV lidar, the system comprising: The data acquisition and time synchronization module is used to acquire laser ranging values, scanning angle, angular velocity, acceleration, position coordinates and pulse count values, and to perform time alignment and marking on the laser ranging values, scanning angle, angular velocity, acceleration and position coordinates based on the pulse count values ​​generated by the field programmable gate array, so as to generate a time-stamped original measurement dataset. The multi-source data interpolation processing module is used to retrieve inertial navigation and positioning data frames that match the laser ranging value from the original time-stamped measurement dataset based on the principle of time proximity, and to use a linear interpolation algorithm to calculate the instantaneous attitude and instantaneous three-dimensional spatial coordinates at each moment. The three-dimensional point cloud coordinate calculation module is used to calculate the three-dimensional coordinates of the sensor coordinate system based on the laser ranging value and the scanning angle, and perform multi-level coordinate transformation in combination with the installation error matrix, the instantaneous attitude and the instantaneous three-dimensional spatial coordinates to generate a three-dimensional laser point cloud dataset in the world geodetic coordinate system. The geological body filtering and modeling module is used to map the three-dimensional laser point cloud dataset to a grid to extract potential ground seed points, filter ground points by calculating local slope values ​​and combining them with slope thresholds, and construct a three-dimensional digital model of the geological body containing irregular triangular meshes and slope and aspect information using the Delaunay triangulation algorithm.

[0015] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, a field-programmable gate array (FPGA) is used to uniformly time-label multi-source sensor data. Adjacent time frame data is retrieved and instantaneous attitude and position are calculated using a linear interpolation algorithm. Based on the installation error matrix, the laser ranging values ​​are accurately converted into a three-dimensional laser point cloud to replace manual contact measurement. The point cloud is mapped to a two-dimensional grid to extract potential ground seed points. The ratio of neighborhood elevation difference to horizontal distance is calculated to generate local slope. Non-ground noise points are filtered and eliminated based on the slope threshold to lock the real ground points. A geological model is constructed using the Delaunay triangulation algorithm. High-density surface data can be obtained without personnel entering dangerous areas, effectively solving the problem of lost details in complex terrain and significantly improving the realism of geological modeling. Attached Figure Description

[0016] Figure 1 This is the main flowchart of the geological body modeling method based on UAV lidar of the present invention; Figure 2 This is a flowchart of the original measurement data acquisition and preprocessing process of this invention; Figure 3 This is a flowchart of the instantaneous attitude and coordinate interpolation calculation of the present invention; Figure 4 This is a flowchart of the three-dimensional laser point cloud generation and coordinate transformation process of the present invention; Figure 5 This is a flowchart of the ground point extraction and geological body modeling process of the present invention. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this invention clearer, the software-based technical solution is described in detail below with reference to system architecture diagrams and embodiments. It should be understood that the specific embodiments described herein are only for explaining the technical solutions of this invention and do not constitute a limitation on the scope of protection.

[0018] In the description of this invention, the system architecture relationships or data processing flows indicated by terms such as "layer," "module," "interface," "data flow," "client," and "server" are all defined based on the architecture diagram or flowchart corresponding to the embodiments. This way of describing is only used to clearly illustrate the logical relationships between the elements in the technical solution, and not to limit the physical deployment form. The term "multiple" includes two or more technical units, including but not limited to multiple data nodes, processing threads, service instances, or functional components and other scalable elements. The specific number is determined according to the actual business scenario and needs to be specifically specified.

[0019] Please see Figure 1 and Figure 2 This invention provides a technical solution: a geological body modeling method based on UAV lidar, comprising the following steps: S1: Acquire laser ranging value, scanning angle, angular velocity, acceleration, position coordinates and pulse count value; read the pulse count value generated by the field programmable gate array; perform time-stamped marking processing on the laser ranging value, scanning angle, angular velocity, acceleration and position coordinates; generate a time-stamped original measurement dataset. The timestamped raw measurement dataset includes laser ranging values, scan angles, angular velocities, accelerations, and position coordinates, all marked with pulse counts at the same time.

[0020] In a specific implementation scenario where a UAV equipped with a lidar system performs geological surveying tasks, after the system powers on, the Field Programmable Gate Array (FPGA) module immediately initializes its internal counter, which accumulates at a clock frequency of 100MHz. When the lidar sensor's transmitting module emits a laser pulse, it simultaneously sends a hardware interrupt signal to the FPGA. Upon detecting the rising edge of this interrupt signal, the FPGA immediately latches the current pulse count value and writes it into the register of the high-speed FIFO buffer. Simultaneously, the onboard inertial navigation system (INS) outputs angular velocity, acceleration, and position coordinate data at a frequency of 200Hz, and the Global Navigation Satellite System (GNSS) receiver outputs a pulse-per-second (PPS) signal and corresponding Coordinated Universal Time (UTC) information at a frequency of 1Hz. The main control computer reads the laser ranging value and scan angle data via a serial port and reads the latched pulse count value from the FPGA.

[0021] During time synchronization and dataset generation operations, the system establishes a linear mapping between the FPGA pulse count value at the arrival time of the PPS signal and the corresponding UTC time. Assuming at time... The PPS signal is received, and the FPGA count value is [value missing]. The FPGA count value read at the current laser emission moment is Then the precise absolute time of laser emission Through calculation formula It is concluded that, among them This is the FPGA clock frequency (100MHz). The system will calculate... As a timestamp, it is written into the header of the raw measurement data packet. For inertial navigation data, the system also records the absolute timestamp of its acquisition time. Next, the system allocates a contiguous storage space in memory and constructs an array of structures. Each structure element contains: a 64-bit double-precision floating-point timestamp, a 32-bit floating-point laser ranging value (unit: meters), a 32-bit floating-point scan angle (unit: degrees), a 32-bit floating-point angular velocity triaxial component (unit: radians / second), a 32-bit floating-point acceleration triaxial component (unit: meters / second squared), and 64-bit double-precision floating-point position coordinates (longitude, latitude, and elevation). As the data stream continues to input, the system merges and aligns these heterogeneous data according to the timestamp order to form a timestamped raw measurement dataset, which is then stored in binary file format on the onboard solid-state drive, completing the data acquisition process.

[0022] Please see Figure 1 and Figure 3 S2: Retrieve angular velocity, acceleration, and position coordinate data frames that are adjacent to the time of the laser ranging value from the original time-stamped measurement dataset, and use the data frames to perform a linear interpolation algorithm to generate instantaneous attitude and instantaneous three-dimensional spatial coordinates; The instantaneous attitude is a quaternion of the body attitude calculated based on linear interpolation; Instantaneous three-dimensional spatial coordinates include interpolated longitude, latitude, and elevation values; The process of generating instantaneous attitude and instantaneous three-dimensional spatial coordinates specifically includes: Read the laser emission time timestamp corresponding to each laser ranging value in the original measurement dataset with timestamps; Using the laser emission time stamp as the retrieval benchmark, a bidirectional search is performed in the time series of the original measurement dataset with timestamps to locate the preceding inertial navigation record frame at the time before the laser emission time stamp and the following inertial navigation record frame at the time after the laser emission time stamp. Extract the position coordinates, angular velocity, and acceleration contained in the preceding and subsequent inertial navigation recording frames, respectively, and calculate the time difference ratio between the laser emission time timestamp and the preceding inertial navigation recording frame. Based on the time difference ratio, linear position weighting calculation is performed on the position coordinates of the preceding and subsequent inertial navigation recording frames to generate instantaneous three-dimensional spatial coordinates; The attitude angle of the aircraft is calculated by fusing angular velocity and acceleration. Spherical linear interpolation is performed on the attitude quaternions of adjacent moments using the time difference ratio to generate the instantaneous attitude corresponding to the timestamp of the laser emission moment.

[0023] In the data processing phase, the computer first loads the timestamped raw measurement dataset stored on the solid-state drive into memory. The program then iterates through each lidar measurement record in the dataset, extracting its associated laser emission timestamp. Using this timestamp as the retrieval reference, a binary search algorithm is performed in the time series index table of the inertial navigation data stream. The system first locates the time index that is less than and closest to... Pre-inertial navigation recording frame Read its timestamp Subsequently, the time index was located to be greater than and closest to Post-inertial navigation recording frames Read its timestamp .

[0024] At this point, the system extracts... Position coordinates in With attitude quaternions ,as well as Position coordinates in With attitude quaternions To ensure calculation accuracy, the position coordinates here are in metric coordinates under the Universal Transverse Mercator (UTM) coordinate system. Assuming... , (Corresponding to east (X), north (Y), and sky (Z) coordinates respectively, unit: meters), attitude quaternion (Unit quaternion, indicating no rotation); (i.e., 200Hz sampling interval) Posture Quaternion (Indicates approximately 5 degrees of rotation around the X-axis); Laser emission time .

[0025] The system first calculates the time difference ratio. The calculation process is as follows .

[0026] Based on this ratio The system performs linear position weighted calculations on the three components of the three-dimensional position coordinates to generate instantaneous three-dimensional spatial coordinates. For the X-axis component (eastward coordinates), the calculation formula is as follows: Substituting the numerical values, we get Meters. For the Y-axis component (northward coordinates), the calculation formula is as follows: Substituting the numerical values, we get Meters. For the Z-axis component (elevation coordinates), the calculation formula is: Substituting the numerical values, we get Meters. The final instantaneous three-dimensional spatial coordinates are: .

[0027] For attitude data processing, the system utilizes time difference ratios. attitude quaternions at adjacent time points and Perform spherical linear interpolation (Slerp). First, calculate the dot product of the two quaternions. Calculate the included angle Since the dot product is greater than 0, shortest path interpolation is used directly. In the interpolation calculation, the weighting factors are as follows: and Substitute Calculated , Instantaneous attitude quaternion According to calculations, The result indicates that in At that moment, the machine underwent a roll rotation of approximately 2 degrees relative to the world coordinate system.

[0028] Please see Figure 1 and Figure 4 S3: Calculate the three-dimensional coordinates of the sensor coordinate system based on the laser ranging value and scanning angle. Based on the installation error matrix, instantaneous attitude and instantaneous three-dimensional spatial coordinates, transform the three-dimensional coordinates of the sensor coordinate system to generate a three-dimensional laser point cloud dataset. The 3D laser point cloud dataset consists of point cloud data in the world geodetic coordinate system after installation error compensation and coordinate system transformation. The calculation process of the three-dimensional coordinates of the sensor coordinate system specifically includes: Obtain the laser ranging values ​​and scanning angles recorded in the original measurement dataset with timestamps, and read the internal optical parameters of the lidar device; The optical path of the laser ranging value is corrected based on the internal optical parameters. The corrected laser ranging value and the scanning angle are decomposed into component data of the laser beam in the rectangular coordinate system defined by the sensor itself using trigonometric function relationships. Based on the horizontal, vertical and longitudinal components obtained from the decomposition, coordinate vectors are constructed to generate the three-dimensional coordinates of the sensor coordinate system. The calculation of the three-dimensional coordinates of the sensor coordinate system follows the formula below: ; in, This represents the generated three-dimensional coordinate vector of the sensor coordinate system. Represents the horizontal axis coordinate value in the sensor coordinate system. This represents the vertical axis coordinate value in the sensor coordinate system. Represents the vertical axis coordinate value in the sensor coordinate system. Represents the laser ranging value. Represents the scanning angle. This represents the fixed pitch deviation angle of the laser emitter relative to the sensor reference plane; The process of generating a 3D laser point cloud dataset specifically includes: Obtain the pre-calibrated installation error matrix, which includes the rotation and translation parameters of the lidar relative to the center of the UAV body; Read the instantaneous attitude and convert it into a rotation matrix of the body coordinate system relative to the geodetic coordinate system; read the instantaneous three-dimensional spatial coordinates as the translation vector of the body center in the geodetic coordinate system. Multiply the three-dimensional coordinates of the sensor coordinate system by the installation error matrix to transform them into intermediate transition coordinates in the body coordinate system. Multiply the intermediate transition coordinates by the rotation matrix and add the translation vector to complete the mapping transformation from the sensor's local space to the world's ground space; The calculation results corresponding to all laser ranging values ​​are iterated through, and all the transformed coordinate points are stacked and stored in chronological order to generate a 3D laser point cloud dataset.

[0029] In the process of generating a 3D laser point cloud dataset, the system first obtains a pre-calibrated installation error matrix. With translation parameters As shown in Table 1, this table records the installation deviation of the lidar sensor relative to the center of the UAV body.

[0030] Table 1. Calibration Parameters for LiDAR Installation Error; Referring to Table 1, the system reads the translation vector. Since the rotational deviation angles in Table 1 are extremely small (all less than 0.2 degrees), in constructing the installation error matrix... When, it is approximated as an identity matrix Further explanation will follow.

[0031] Next, the system performs coordinate calculations for each laser point. First, it reads the laser ranging values ​​recorded in the timestamped raw measurement dataset. With scanning angle Assuming the current time... rice, (Approximately 0.5236 radians). Read the fixed pitch deviation angle of the laser emitter relative to the sensor reference plane. , set as After correcting the optical path of the laser ranging value based on internal optical parameters, the system uses trigonometric relationships to decompose the corrected laser ranging value and the scanning angle, generating a three-dimensional coordinate system for the sensor. .

[0032] The calculation of the three-dimensional coordinates of the sensor coordinate system follows the formula below: ; in, This represents the generated three-dimensional coordinate vector of the sensor coordinate system. Represents the horizontal axis coordinate value in the sensor coordinate system. This represents the vertical axis coordinate value in the sensor coordinate system. Represents the vertical axis coordinate value in the sensor coordinate system. Represents the laser ranging value. Represents the scanning angle. This represents the fixed pitch deviation angle of the laser emitter relative to the sensor reference plane.

[0033] Substitute the above example values ​​into the formula to calculate: rice; rice; Rice. That is... .

[0034] Subsequently, the system performs a coordinate system transformation. First, the three-dimensional coordinates of the sensor coordinate system are multiplied by the installation error matrix. And add translation vector The intermediate transition coordinates in the body coordinate system are obtained. The calculation process is as follows: .

[0035] Next, the system multiplies the intermediate transition coordinates on the left by the instantaneous attitude quaternion. The resulting rotation matrix is ​​then added to the instantaneous three-dimensional spatial coordinates (i.e., the translation vector). Referring to the calculation results in step S2, Corresponding to a roll angle of approximately 2 degrees, its rotation matrix Acting on This will cause a rotation transformation of the Y-axis and Z-axis components. The rotated coordinate vector The calculation is as follows: rice; rice; Meters. Finally, add the translation vector. Final world coordinates The system iterates through all the calculation results corresponding to the laser ranging values, stacks and stores all the transformed coordinate points in chronological order, and generates a 3D laser point cloud dataset.

[0036] Please see Figure 1 and Figure 5 S4: Map the 3D laser point cloud dataset to a 2D grid to extract potential ground seed points, calculate the ratio of the elevation difference and horizontal distance between the potential ground seed points and their neighbors to generate local slope values, filter local slope values ​​based on slope thresholds to determine ground points, and input the ground points into the Delaunay triangulation algorithm to generate a 3D digital model of the geological body. The process of extracting potential ground seed points specifically includes: Read the extreme values ​​of the planar coordinates of all point cloud data in the 3D laser point cloud dataset to determine the horizontal rectangular boundary range of the point cloud coverage. Based on the preset grid resolution parameters, the horizontal rectangular boundary range is divided into non-overlapping regular two-dimensional grid matrices; Traverse each laser point in the 3D laser point cloud dataset and map the laser point index to the corresponding 2D grid cell based on the laser point's planar coordinates. For each two-dimensional grid cell containing laser points, retrieve the elevation values ​​of all laser points within that two-dimensional grid cell, and filter out the laser point with the minimum elevation value through comparison calculation; The laser point with the lowest elevation value is marked as the representative low point of the two-dimensional grid cell, and the representative low points of all two-dimensional grid cells are aggregated to generate potential ground seed points. The process of generating local slope values ​​specifically includes: Construct a KD-tree spatial index structure for potential ground seed points. For each potential ground seed point, use the KD-tree spatial index structure to search for its nearest neighbor set within a set radius. Traverse each neighboring point in the nearest neighbor set and extract the elevation and planar coordinate components of the potential ground seed point and that neighboring point respectively. Calculate the absolute value of the elevation difference between the potential ground seed point and its neighboring point, and calculate the Euclidean distance between the potential ground seed point and its neighboring point projected onto the horizontal plane. The slope tangent is calculated based on the ratio of the absolute value of the elevation difference to the Euclidean distance. The maximum value among the calculation results of all neighboring points is selected to generate the local slope value. The calculation of local slope values ​​follows the formula below: ; in, This represents the generated local slope value. The set of indices representing the nearest neighbor set. The elevation values ​​representing potential ground seed points. Representing the Elevation values ​​of neighboring points, The lateral plane coordinates of potential ground seed points. The vertical plane coordinates of the potential ground seed points are represented. Representing the The horizontal plane coordinates of the neighboring points Representing the The vertical plane coordinate values ​​of the neighboring points; The process of determining ground points specifically includes: Obtain the preset slope threshold, which is set based on the average roughness of the geological type of the survey area; Read the local slope value corresponding to each potential ground seed point, and compare the local slope value with the slope threshold one by one; If the local slope value is less than or equal to the slope threshold, it is determined that the local terrain of the potential ground seed point is gentle, and the attribute identifier of the potential ground seed point is modified to a confirmed ground attribute. If the local slope value is greater than the slope threshold, the potential ground seed point is determined to be a vegetation or non-ground obstacle interference point, and the potential ground seed point is removed from the candidate set. Extract all point data marked as having certain ground attributes, perform outlier denoising on the point data with certain ground attributes, and generate ground points; The process of generating a three-dimensional digital model of a geological body specifically includes: Project the ground points onto a two-dimensional plane to construct the convex hull boundary of the point set; The Deloitte triangulation algorithm is called to connect the ground points as vertices to generate a non-overlapping triangular topology network that covers the entire convex hull boundary, thus generating an irregular triangular network. Traverse each triangular facet in the irregular triangular mesh, extract the three-dimensional coordinates of the three vertices, and calculate the normal vector of the triangular facet through vector cross product operation; The slope parameters are calculated based on the vertical component of the normal vector, and the aspect parameters are calculated based on the horizontal projection of the normal vector. The irregular triangular mesh, triangular patch normal vectors, geological body slope parameters and aspect parameters are encapsulated to generate a three-dimensional digital model of the geological body. The three-dimensional digital model of the geological body includes an irregular triangular mesh generated based on ground points, triangular patch normal vectors, geological body slope parameters, and slope aspect parameters.

[0037] In the stage of generating a 3D digital model of the geological body, the system first reads the extreme values ​​of the planar coordinates of all point cloud data in the 3D laser point cloud dataset to determine the horizontal rectangular boundary range covered by the point cloud. Based on a preset grid resolution parameter (e.g., 1.0 meter), the horizontal rectangular boundary range is divided into a non-overlapping regular 2D grid matrix. Each laser point in the 3D laser point cloud dataset is traversed, and the laser point index is mapped to the corresponding 2D grid cell based on its planar coordinates. For each 2D grid cell containing a laser point, the elevation values ​​of all laser points within that grid cell are retrieved, and the laser point with the lowest elevation value is selected through comparison and marked as a potential ground seed point.

[0038] Subsequently, the system constructs a KD-tree spatial index structure for potential ground seed points. For each potential ground seed point, the KD-tree spatial index structure is used to search for its nearest neighbor set within a defined radius (e.g., 1.5 meters). The coordinates calculated in step S3 are then used as the basis for this search. The current seed point The system found two points in the neighborhood. and The specific values ​​are shown in Table 2.

[0039] Table 2 Spatial data of seed points and their neighboring points; Referring to Table 2, the system calculates the ratio of the absolute value of the elevation difference between the potential ground seed point and the neighboring point to the Euclidean distance, and generates local slope values.

[0040] The calculation of local slope values ​​follows the formula below: ; in, This represents the generated local slope value. The set of indices representing the nearest neighbor set. The elevation values ​​representing potential ground seed points. Representing the Elevation values ​​of neighboring points, The lateral plane coordinates of potential ground seed points. The vertical plane coordinates of the potential ground seed points are represented. Representing the The horizontal plane coordinates of the neighboring points Representing the The vertical plane coordinate values ​​of the neighboring points.

[0041] Substitute the data from Table 2 into a practical example calculation: for neighborhood points Horizontal distance meters; elevation difference meters; slope value For neighboring points Horizontal distance meters; elevation difference meters; slope value The system selects the maximum value from the calculation results to generate a local slope value. .

[0042] The process of determining ground points specifically includes: obtaining a preset slope threshold, which is set based on the average roughness of the geological type of the survey area; for example, through statistical analysis of flat areas, the slope threshold is set to 0.80. The calculated local slope values ​​are then read. This was compared to a slope threshold of 0.80. Because... If a potential ground seed point is determined to be a vegetation or non-ground obstacle interference point (such as the top of a bush), it will be removed from the candidate set. If the local slope value of a seed point is less than or equal to 0.80, the attribute identifier of the potential ground seed point will be modified to a confirmed ground attribute.

[0043] Extract all point data marked as having certain ground properties, perform outlier denoising on these points, and generate ground points. The specific process of generating the 3D digital model of the geological body includes: projecting ground points onto a 2D plane to construct the convex hull boundary of the point set; using the Delaunay triangulation algorithm, connecting ground points as vertices to generate a non-overlapping triangular topology network that covers the entire convex hull boundary, thus generating an irregular triangular network; traversing each triangular facet in the irregular triangular network, extracting the 3D coordinates of the three vertices, and calculating the normal vector of the facet using vector cross product; calculating the slope parameter based on the vertical component of the normal vector, and calculating the aspect parameter based on the horizontal projection direction of the normal vector; encapsulating the irregular triangular network, triangular facet normal vectors, geological body slope parameters, and aspect parameters to generate the 3D digital model of the geological body.

[0044] A geological body modeling system based on UAV LiDAR, the system being used to execute the aforementioned geological body modeling method based on UAV LiDAR, the system comprising: The data acquisition and time synchronization module is used to acquire laser ranging values, scanning angle, angular velocity, acceleration, position coordinates and pulse count values, and to perform time alignment and marking on the laser ranging values, scanning angle, angular velocity, acceleration and position coordinates based on the pulse count values ​​generated by the field programmable gate array, generating a time-stamped original measurement dataset; The multi-source data interpolation processing module is used to retrieve inertial navigation and positioning data frames that match the laser ranging values ​​from the original time-stamped measurement dataset based on the principle of time proximity, and to use a linear interpolation algorithm to solve the instantaneous attitude and instantaneous three-dimensional spatial coordinates at each moment. The 3D point cloud coordinate calculation module is used to calculate the 3D coordinates of the sensor coordinate system based on the laser ranging value and the scanning angle, and to perform multi-level coordinate transformations in combination with the installation error matrix, instantaneous attitude and instantaneous 3D spatial coordinates to generate a 3D laser point cloud dataset in the world geodetic coordinate system. The geological body filtering and modeling module is used to map the 3D laser point cloud dataset to a grid to extract potential ground seed points. It filters ground points by calculating local slope values ​​and combining them with slope thresholds, and uses the DeLorean triangulation algorithm to construct a 3D digital model of the geological body containing irregular triangular meshes and slope and aspect information.

[0045] The above embodiments illustrate preferred embodiments of the present invention. Any equivalent adjustments to the technical solution based on software engineering methods are within the scope of protection, including but not limited to: implementing algorithm logic using different programming languages, refactoring functional modules into services, adjusting data interaction protocols, and optimizing resource scheduling strategies. Any implementation scheme derived from reasonable modifications to the data processing flow, service call chain, or system architecture layer without departing from the core technology of the present invention should be considered within the scope of protection defined by the claims of the present invention.

Claims

1. A geological body modeling method based on UAV lidar, characterized in that, Includes the following steps: S1: Acquire laser ranging value, scanning angle, angular velocity, acceleration, position coordinates and pulse count value; read the pulse count value generated by field programmable gate array; perform time-stamping processing on the laser ranging value, the scanning angle, the angular velocity, the acceleration and the position coordinates to generate a time-stamped original measurement dataset; S2: Retrieve the angular velocity, acceleration, and position coordinate data frames that are adjacent to the time of the laser ranging value from the original time-stamped measurement dataset, and use the data frames to perform a linear interpolation algorithm to generate instantaneous attitude and instantaneous three-dimensional spatial coordinates; S3: Calculate the three-dimensional coordinates of the sensor coordinate system based on the laser ranging value and the scanning angle; and generate a three-dimensional laser point cloud dataset by transforming the three-dimensional coordinates of the sensor coordinate system based on the installation error matrix, the instantaneous attitude, and the instantaneous three-dimensional spatial coordinates. S4: Map the three-dimensional laser point cloud dataset to a two-dimensional grid to extract potential ground seed points, calculate the ratio of the elevation difference and horizontal distance between the potential ground seed points and their neighbors to generate local slope values, filter the local slope values ​​based on slope thresholds to determine ground points, and input the ground points into the Delaunay triangulation algorithm to generate a three-dimensional digital model of the geological body.

2. The geological body modeling method based on UAV lidar according to claim 1, characterized in that, The timestamped raw measurement dataset includes the laser ranging value, the scanning angle, the angular velocity, the acceleration, and the position coordinates, all marked with pulse count values ​​at the same time. The instantaneous attitude is a quaternion of the body attitude calculated based on linear interpolation; The instantaneous three-dimensional spatial coordinates include interpolated longitude, latitude, and elevation values; The three-dimensional laser point cloud dataset consists of point cloud data in the world geodetic coordinate system after installation error compensation and coordinate system transformation. The three-dimensional digital model of the geological body includes an irregular triangular mesh generated based on the ground points, triangular patch normal vectors, geological body slope parameters, and slope aspect parameters.

3. The geological body modeling method based on UAV lidar according to claim 2, characterized in that, The process of generating the instantaneous attitude and the instantaneous three-dimensional spatial coordinates specifically includes: Read the laser emission time timestamp corresponding to each laser ranging value in the timestamped raw measurement dataset; Using the laser emission time timestamp as the retrieval reference, a bidirectional search is performed in the time series of the original measurement dataset with timestamps to locate the preceding inertial navigation recording frame at the time before the laser emission time timestamp and the following inertial navigation recording frame at the time after the laser emission time timestamp. Extract the position coordinates, angular velocity and acceleration contained in the preceding inertial navigation recording frame and the following inertial navigation recording frame respectively, and calculate the time difference ratio of the laser emission time timestamp relative to the preceding inertial navigation recording frame; Based on the time difference ratio, linear position weighting calculation is performed on the position coordinates of the preceding inertial navigation recording frame and the following inertial navigation recording frame to generate the instantaneous three-dimensional spatial coordinates; The body attitude angle is calculated by fusing the angular velocity and the acceleration, and spherical linear interpolation is performed on the attitude quaternions of adjacent moments using the time difference ratio to generate the instantaneous attitude corresponding to the timestamp of the laser emission moment.

4. The geological body modeling method based on UAV lidar according to claim 3, characterized in that, The calculation process of the three-dimensional coordinates of the sensor coordinate system specifically includes: Obtain the laser ranging value and the scanning angle recorded in the original measurement dataset with timestamps, and read the internal optical parameters of the lidar device; The laser ranging value is optically corrected based on the internal optical parameters. The corrected laser ranging value and the scanning angle are decomposed into component data of the laser beam in the rectangular coordinate system defined by the sensor itself using trigonometric function relationships. Based on the horizontal, vertical, and longitudinal components obtained from the decomposition, coordinate vectors are constructed to generate the three-dimensional coordinates of the sensor coordinate system. The calculation of the three-dimensional coordinates of the sensor coordinate system follows the formula below: ； in, This represents the generated three-dimensional coordinate vector of the sensor coordinate system. Represents the horizontal axis coordinate value in the sensor coordinate system. This represents the vertical axis coordinate value in the sensor coordinate system. Represents the vertical axis coordinate value in the sensor coordinate system. This represents the laser ranging value. Represents the scanning angle, This represents the fixed pitch deviation angle of the laser emitter relative to the sensor reference plane.

5. The geological body modeling method based on UAV lidar according to claim 4, characterized in that, The generation process of the three-dimensional laser point cloud dataset specifically includes: Obtain the pre-calibrated installation error matrix, which includes the rotation and translation parameters of the lidar relative to the center of the UAV body; The instantaneous attitude is read and converted into a rotation matrix of the body coordinate system relative to the geodetic coordinate system. The instantaneous three-dimensional spatial coordinates are read as the translation vector of the body center in the geodetic coordinate system. Multiply the three-dimensional coordinates of the sensor coordinate system by the installation error matrix to transform them into intermediate transition coordinates in the body coordinate system. Multiply the intermediate transition coordinates by the rotation matrix and add the translation vector to complete the mapping transformation from the sensor's local space to the world's ground space; The calculation results corresponding to all the laser ranging values ​​are traversed, and all the transformed coordinate points are stacked and stored in chronological order to generate the three-dimensional laser point cloud dataset.

6. The geological body modeling method based on UAV lidar according to claim 5, characterized in that, The extraction process of the potential ground seed points specifically includes: Read the extreme values ​​of the planar coordinates of all point cloud data in the three-dimensional laser point cloud dataset to determine the horizontal rectangular boundary range of the point cloud coverage. Based on preset grid resolution parameters, the horizontal rectangular boundary range is divided into non-overlapping regular two-dimensional grid matrices; Traverse each laser point in the three-dimensional laser point cloud dataset, and map the laser point index to the corresponding two-dimensional grid cell according to the planar coordinates of the laser point; For each two-dimensional grid cell containing the laser point, retrieve the elevation values ​​of all laser points within that two-dimensional grid cell, and filter out the laser point with the minimum elevation value through comparison calculation; The laser point with the minimum elevation value is marked as the representative low point of the two-dimensional grid cell, and the representative low points of all the two-dimensional grid cells are aggregated to generate the potential ground seed point.

7. The geological body modeling method based on UAV lidar according to claim 6, characterized in that, The process of generating the local slope value specifically includes: Construct a KD tree spatial index structure for the potential ground seed points, and for each potential ground seed point, use the KD tree spatial index structure to search for its nearest neighbor set within a set radius. Traverse each neighboring point in the nearest neighbor set and extract the elevation and planar coordinate components of the potential ground seed point and the neighboring point respectively. Calculate the absolute value of the elevation difference between the potential ground seed point and its neighboring point, and calculate the Euclidean distance between the potential ground seed point and its neighboring point projected onto the horizontal plane. The slope tangent is calculated based on the ratio of the absolute value of the elevation difference to the Euclidean distance. The maximum value among the calculation results of all the neighboring points is selected to generate the local slope value. The calculation of the local slope value follows the following formula: ； in, This represents the generated local slope value. The set of indices representing the nearest neighbor set, This represents the elevation value of the potential ground seed point. Representing the The elevation values ​​of the neighboring points. The lateral plane coordinates of the potential ground seed points are represented. The longitudinal plane coordinates of the potential ground seed points are represented. Representing the The lateral plane coordinates of the neighboring points. Representing the The vertical plane coordinate values ​​of the neighboring points.

8. The geological body modeling method based on UAV lidar according to claim 7, characterized in that, The process of determining the ground points specifically includes: Obtain the preset slope threshold, which is set based on the average roughness of the geological type of the survey area; Read the local slope value corresponding to each potential ground seed point, and compare the local slope value with the slope threshold value one by one; If the local slope value is less than or equal to the slope threshold, it is determined that the local terrain of the potential ground seed point is flat, and the attribute identifier of the potential ground seed point is modified to a confirmed ground attribute. If the local slope value is greater than the slope threshold, the potential ground seed point is determined to be a vegetation or non-ground obstacle interference point, and the potential ground seed point is removed from the candidate set. Extract all point data marked with the confirmed ground attribute, perform outlier denoising processing on the point data of the confirmed ground attribute, and generate the ground point.

9. The geological body modeling method based on UAV lidar according to claim 8, characterized in that, The process of generating the three-dimensional digital model of the geological body specifically includes: Project the ground points onto a two-dimensional plane to construct the convex hull boundary of the point set; The Deloitte triangulation algorithm is called to connect the ground points as vertices to generate a non-overlapping triangular topology network that covers the entire convex hull boundary, thus generating an irregular triangular network. Traverse each triangular facet in the irregular triangular network, extract the three-dimensional coordinates of the three vertices, and calculate the normal vector of the triangular facet through vector cross product operation; The slope parameter is calculated based on the vertical component of the normal vector, and the aspect parameter is calculated based on the horizontal projection direction of the normal vector. The irregular triangular mesh, the normal vector of the triangular facet, the slope parameters of the geological body, and the aspect parameters are encapsulated to generate a three-dimensional digital model of the geological body.

10. A geological body modeling system based on UAV lidar, characterized in that, The system is used to implement the geological body modeling method based on UAV lidar as described in any one of claims 1-9, and the system includes: The data acquisition and time synchronization module is used to acquire laser ranging values, scanning angle, angular velocity, acceleration, position coordinates and pulse count values, and to perform time alignment and marking on the laser ranging values, scanning angle, angular velocity, acceleration and position coordinates based on the pulse count values ​​generated by the field programmable gate array, so as to generate a time-stamped original measurement dataset. The multi-source data interpolation processing module is used to retrieve inertial navigation and positioning data frames that match the laser ranging value from the original time-stamped measurement dataset based on the principle of time proximity, and to use a linear interpolation algorithm to calculate the instantaneous attitude and instantaneous three-dimensional spatial coordinates at each moment. The three-dimensional point cloud coordinate calculation module is used to calculate the three-dimensional coordinates of the sensor coordinate system based on the laser ranging value and the scanning angle, and perform multi-level coordinate transformation in combination with the installation error matrix, the instantaneous attitude and the instantaneous three-dimensional spatial coordinates to generate a three-dimensional laser point cloud dataset in the world geodetic coordinate system. The geological body filtering and modeling module is used to map the three-dimensional laser point cloud dataset to a grid to extract potential ground seed points, filter ground points by calculating local slope values ​​and combining them with slope thresholds, and construct a three-dimensional digital model of the geological body containing irregular triangular meshes and slope and aspect information using the Delaunay triangulation algorithm.

Images