Unmanned aerial vehicle remote sensing image assisted terrain elevation data inversion calculation method

Abstract

The present application relates to the technical field of remote sensing topographic mapping, in particular to a method for calculating terrain elevation data inversion assisted by unmanned aerial vehicle remote sensing images, comprising: collecting original remote sensing image sequences and laser point cloud data sets of an optical camera and a laser radar carried by an unmanned aerial vehicle, and generating an initial sparse point cloud three-dimensional coordinate set through stereo image matching; after unifying the spatial coordinate systems of the two types of data, extracting the surface reflection intensity and texture gradient features of the ground objects to construct a joint feature matrix, and inputting the joint feature matrix into a pre-trained deep convolutional neural network to output a terrain category probability distribution map; filtering the point cloud by region based on the distribution map to separate the ground and non-ground point clouds, generating a first edition digital elevation model through triangulation encryption interpolation of the ground point cloud, calculating an elevation residual distribution field, and correcting the elevation values by grid. This method fuses multi-source data and multi-dimensional features, accurately distinguishes terrain categories, corrects the numerical deviation of the elevation model, adapts to complex terrain environments, and optimizes the data fitting degree of terrain elevation inversion.

Description

Technical Field

[0001] This invention relates to the field of remote sensing topographic mapping technology, and in particular to a method for inverting and calculating topographic elevation data assisted by UAV remote sensing images. Background Technology

[0002] In conventional UAV terrain elevation data inversion operations, a single data source, either optical remote sensing imagery or lidar point cloud, is typically used for data processing. Digital elevation models are constructed based on traditional fixed threshold filtering and basic geometric interpolation methods, which is currently the mainstream approach for terrain elevation mapping. However, the surface feature information provided by a single data source is relatively limited, only acquiring surface information in a single dimension. It is difficult to simultaneously consider surface texture details and lidar echo reflection information, and its ability to classify complex terrain areas is limited.

[0003] Traditional point cloud filtering methods lack intelligent classification criteria, relying solely on fixed parameters to separate ground and non-ground point clouds, resulting in poor point cloud separation performance in complex terrain scenarios. Raster elevation models generated by interpolating separated point clouds are prone to numerical deviations from actual terrain coordinates. Conventional modeling processes lack a correction step for these elevation deviations, leading to significant discrepancies between the resulting elevation grid values ​​and the actual terrain undulations, failing to meet the application requirements of high-precision terrain elevation inversion. Therefore, it is necessary to integrate two data sources—optical imaging and LiDAR—to mine multi-dimensional surface features for refined terrain classification, while simultaneously using point elevation differences to calibrate and optimize the elevation model grid values. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of existing technologies and propose a method for inverting and calculating terrain elevation data assisted by UAV remote sensing images.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a method for inverting and calculating terrain elevation data assisted by UAV remote sensing images, comprising:

[0006] The original remote sensing image sequence and original laser point cloud data set captured by the optical camera and lidar on the UAV in the target survey area are collected. The original remote sensing image sequence is subjected to stereo image pair matching and an initial sparse point cloud three-dimensional coordinate set is generated.

[0007] The initial sparse point cloud 3D coordinate set is transformed with the original laser point cloud data set into a spatial coordinate system. The surface reflection intensity features and texture gradient features of the ground objects are extracted from the data after the coordinate system transformation, and a joint feature matrix is ​​constructed.

[0008] The joint feature matrix is ​​input into a pre-trained deep convolutional neural network and a preliminary terrain category probability distribution map is output.

[0009] Based on the preliminary terrain category probability distribution map, the original laser point cloud data set is filtered by region to separate ground point clouds from non-ground point clouds;

[0010] The separated ground point cloud was subjected to irregular triangular mesh encryption interpolation to generate the first version of the raster digital elevation model;

[0011] Calculate the elevation residual distribution field between the first version of the raster-format digital elevation model and the initial sparse point cloud three-dimensional coordinate set;

[0012] The first raster digital elevation model is corrected grid by grid using the elevation residual distribution field to generate a second raster digital elevation model. The second raster digital elevation model is then output as the final terrain elevation data inversion calculation result.

[0013] As a further aspect of the present invention, the step of collecting the original remote sensing image sequence and the original laser point cloud data set captured by the optical camera and lidar mounted on the UAV on the target survey area specifically includes:

[0014] Control the drone to perform grid-based full-coverage flight over the target survey area according to the preset route and preset flight altitude;

[0015] During the flight of the UAV, the optical camera is synchronously triggered to acquire ground optical images at fixed time intervals and form the original remote sensing image sequence.

[0016] During the flight of the UAV, the lidar is synchronously triggered to collect three-dimensional spatial coordinate points at a fixed pulse frequency and form the original laser point cloud data set.

[0017] Record the instantaneous position coordinates and attitude angle data of the UAV each time optical images are acquired;

[0018] Record the instantaneous position coordinates and attitude angle data of the UAV each time laser point cloud is acquired;

[0019] The original remote sensing image sequence and the original laser point cloud data set are stored in time-aligned according to the acquisition timestamp.

[0020] As a further aspect of the present invention, the step of performing stereo pair matching on the original remote sensing image sequence and generating an initial sparse point cloud three-dimensional coordinate set specifically includes:

[0021] Extract overlapping image pairs between adjacent flight strips from the original remote sensing image sequence to form a stereo image pair set;

[0022] Feature point detection is performed on each pair of overlapping images in the stereo image pair set to extract scale-invariant feature points;

[0023] Nearest neighbor matching is performed on scale-invariant feature points detected between two adjacent images to establish correspondences of corresponding points.

[0024] The disparity value of each matching point pair is calculated based on the difference in pixel coordinates corresponding to the same point.

[0025] The ground 3D coordinates corresponding to each matching point are calculated by combining the instantaneous position coordinates and attitude angle data of the UAV with the parallax value using the forward intersection method;

[0026] All the calculated 3D ground coordinates are aggregated into the initial sparse point cloud 3D coordinate set.

[0027] As a further aspect of the present invention, the step of performing a spatial coordinate system transformation between the initial sparse point cloud three-dimensional coordinate set and the original laser point cloud data set specifically includes:

[0028] Extract the latitude, longitude, and elevation coordinates of each point in the initial sparse point cloud 3D coordinate set;

[0029] Extract the local projected coordinate system plane coordinates and elevation coordinates of each point in the original laser point cloud data set;

[0030] Select at least three common control points between the initial sparse point cloud 3D coordinate set and the original laser point cloud data set;

[0031] Calculate the seven-parameter transformation matrix between the geodetic coordinate system and the local projected coordinate system based on the aforementioned common control points;

[0032] Transform all points in the initial sparse point cloud 3D coordinate set to the local projected coordinate system according to the seven-parameter transformation matrix;

[0033] The initial sparse point cloud 3D coordinate set after coordinate system transformation is superimposed and stored with the original laser point cloud data set.

[0034] As a further aspect of the present invention, the step of extracting surface reflection intensity features and texture gradient features from the data after coordinate system transformation and constructing a joint feature matrix specifically includes:

[0035] Extract the echo reflection intensity value of each laser point from the original laser point cloud data set after coordinate system transformation;

[0036] The echo reflection intensity value of each laser point is mapped onto a two-dimensional regular grid according to the planar coordinate position of that point to form a reflection intensity feature map;

[0037] Extract the original remote sensing image pixel value corresponding to each point from the initial sparse point cloud three-dimensional coordinate set after coordinate system transformation;

[0038] Calculate the spatial rate of change of each pixel value in the horizontal and vertical directions as texture gradient features;

[0039] The reflection intensity feature value and the texture gradient feature value at the same grid location are concatenated in the order of grid rows and columns to form a two-dimensional feature vector group;

[0040] The two-dimensional feature vectors of all grids are arranged in spatial order to form a three-dimensional tensor as the joint feature matrix.

[0041] As a further aspect of the present invention, the step of inputting the joint feature matrix into a pre-trained deep convolutional neural network and outputting a preliminary terrain category probability distribution map specifically includes:

[0042] Obtain a deep convolutional neural network model that has been pre-trained on a labeled terrain sample dataset;

[0043] The joint feature matrix is ​​cut into multiple fixed-size feature patches according to the spatial coverage of the original remote sensing image sequence;

[0044] Each feature map block is sequentially input into the deep convolutional neural network model for forward propagation calculation;

[0045] The deep convolutional neural network model outputs a probability vector for each feature map block;

[0046] Each component in the probability vector corresponds to a feature patch belonging to a preset terrain category;

[0047] The output probability vectors of all feature map blocks are reassembled into a multi-channel probability map according to their original spatial positions.

[0048] The terrain category with the highest probability at each spatial location in the multi-channel probability map is marked as the preliminary terrain category at that location, and the preliminary terrain category probability distribution map is formed.

[0049] The training process of the deep convolutional neural network model includes: extracting and enhancing features from the joint feature matrix through a multi-scale feature fusion module.

[0050] As a further aspect of the present invention, the step of performing regional filtering on the original laser point cloud data group based on the preliminary terrain category probability distribution map to separate ground point clouds from non-ground point clouds specifically includes:

[0051] Each laser point in the original laser point cloud data set is located at its corresponding position in the preliminary terrain category probability distribution map according to its planar coordinates;

[0052] Read the terrain category with the highest probability at that location and use it as the terrain category label for that laser point;

[0053] When the terrain category label is bare ground or low vegetation, the laser point is marked as a candidate ground point;

[0054] When the terrain category label is building or tall vegetation, the laser point is marked as a candidate non-ground point;

[0055] Progressive morphological filtering is performed on the laser point cloud marked as candidate ground points to remove outliers;

[0056] Clustering segmentation filtering is performed on the laser point cloud labeled as candidate non-ground points to separate the independent non-ground object point cloud;

[0057] The remaining ground points and non-ground points after filtering are stored in the ground point cloud set and the non-ground point cloud set, respectively.

[0058] As a further aspect of the present invention, the step of performing irregular triangular mesh encryption interpolation on the separated ground point cloud and generating a first version of the raster-format digital elevation model specifically includes:

[0059] Extract the planar coordinates and elevation values ​​of all ground points from the ground point cloud set;

[0060] Construct a Delaunay irregular triangulation using the planar coordinates of all ground points as the vertices of triangles;

[0061] Linear interpolation is performed inside each triangle of the Delaunay irregular triangular network according to the elevation values ​​of the three vertices of the triangle;

[0062] Arrange regular grid points within the target survey area according to the preset grid spacing;

[0063] For each regular grid point, locate the irregular triangular network triangle it belongs to and calculate the interpolated elevation value inside the triangle;

[0064] Organize the elevation values ​​of all regular grid points into a two-dimensional matrix according to the grid row and column order;

[0065] Inverse distance weighted interpolation is performed on the blank grid points in the two-dimensional matrix to fill the blanks.

[0066] The two-dimensional matrix after filling the blanks is used as the first version of the raster digital elevation model.

[0067] As a further aspect of the present invention, the step of calculating the elevation residual distribution field between the first raster-format digital elevation model and the initial sparse point cloud three-dimensional coordinate set specifically includes:

[0068] The elevation value of each regular grid point in the first version of the raster digital elevation model is used as the reference elevation value.

[0069] The elevation value of each sparse point in the initial sparse point cloud three-dimensional coordinate set is used as the observed elevation value.

[0070] For each sparse point, find the grid point in the first version of the raster digital elevation model to which the plane coordinates of that sparse point belong;

[0071] The difference between the observed elevation value of the sparse point and the reference elevation value of the grid point is calculated as the elevation residual of the point.

[0072] Allocate the elevation residuals of all sparse points to the corresponding grid cells according to their planar coordinate positions;

[0073] The arithmetic mean of multiple elevation residuals assigned to the same grid cell is calculated and used as the representative residual value of that grid cell;

[0074] An incomplete residual field is formed by combining grid cells that represent residual values ​​with null grid cells that do not have sparse points.

[0075] Kriging spatial interpolation is performed on the null grid cells in the incomplete residual field to generate a complete gridded residual field as the elevation residual distribution field.

[0076] As a further aspect of the present invention, the step of using the elevation residual distribution field to perform grid-by-grid elevation value correction on the first raster digital elevation model to generate a second raster digital elevation model, and outputting the second raster digital elevation model as the final terrain elevation data inversion calculation result, specifically includes:

[0077] The residual value of each grid cell in the elevation residual distribution field is read from that grid cell;

[0078] The reference elevation value of the same grid cell in the first version of the raster digital elevation model is algebraically summed with the residual value of that grid cell.

[0079] The summation result is written as the corrected elevation value of the grid cell into the corresponding position of a new blank two-dimensional matrix;

[0080] The steps of reading the residual value and calculating the algebraic summation are repeated for each grid cell in the blank two-dimensional matrix until all grid cells have been processed.

[0081] The two-dimensional matrix containing all the corrected elevation values ​​is used as the second version of the raster-format digital elevation model.

[0082] Export the second version of the raster digital elevation model as a file according to the standard geographic information system raster format;

[0083] The exported file is output as the final terrain elevation data inversion calculation result.

[0084] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0085] The initial sparse point cloud 3D coordinate set generated by matching stereo image pairs from the UAV optical camera is transformed into a single spatial coordinate system with the original laser point cloud data set. Simultaneously, surface reflection intensity features and texture gradient features are extracted, and a joint feature matrix is ​​constructed. This matrix is ​​then fed into a pre-trained deep convolutional neural network to output a terrain category probability distribution map. The multi-dimensional surface features complement and fuse each other, broadening the representation range of surface information and overcoming the limitations of single features in terrain recognition. The network model autonomously learns feature relationships, refining the category discrimination of different landform regions. Subsequent point cloud regional filtering can match the actual terrain distribution to perform classification operations.

[0086] The separated ground point cloud is subjected to irregular triangular mesh densification interpolation to generate a raster digital elevation model. The elevation residual distribution field between this model and the initial sparse point cloud 3D coordinate set is solved. The residual distribution field is used to perform grid-by-grid elevation value correction processing on the elevation model. Triangular mesh densification interpolation can refine the terrain transition details between discrete point clouds, and the residual distribution field can intuitively reflect the degree of deviation of elevation values ​​in each region. Grid-by-grid calibration of elevation values ​​smooths out the numerical offset of local grids, weakens the terrain numerical distortion generated during the interpolation modeling process, and makes the rasterized elevation data fit the real undulation of the ground surface. Attached Figure Description

[0087] Figure 1 This is a flowchart of the method for inverting terrain elevation data assisted by UAV remote sensing images as described in this invention;

[0088] Figure 2 A flowchart for data acquisition and time alignment;

[0089] Figure 3 This is a flowchart for filtering and separating ground and non-ground point clouds in different regions. Detailed Implementation

[0090] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0091] See Figure 1This invention provides a method for inverting and calculating terrain elevation data assisted by UAV remote sensing images. The specific method includes:

[0092] The system collects raw remote sensing image sequences and raw laser point cloud data sets from the target survey area captured by the optical camera and lidar mounted on an UAV. Stereo pair matching is performed on the raw remote sensing image sequences to generate an initial sparse point cloud 3D coordinate set. A spatial coordinate system transformation is performed between the initial sparse point cloud 3D coordinate set and the raw laser point cloud data set. Surface reflectance intensity features and texture gradient features are extracted from the transformed data, and a joint feature matrix is ​​constructed. The joint feature matrix is ​​input into a pre-trained deep convolutional neural network, which outputs a preliminary terrain category probability distribution map. Based on... The preliminary terrain category probability distribution map is used to perform regional filtering on the original laser point cloud data set to separate the ground point cloud from the non-ground point cloud; the separated ground point cloud is then subjected to irregular triangular mesh densification interpolation to generate the first version of the raster digital elevation model; the elevation residual distribution field between the first version of the raster digital elevation model and the initial sparse point cloud 3D coordinate set is calculated; the elevation residual distribution field is used to perform grid-by-grid elevation value correction on the first version of the raster digital elevation model to generate the second version of the raster digital elevation model, and the second version of the raster digital elevation model is output as the result of terrain elevation data inversion calculation.

[0093] In one embodiment of the present invention, the process of acquiring the original remote sensing image sequence and the original laser point cloud data set captured by the optical camera and lidar on the target survey area by the UAV is as follows: (See...) Figure 2The system controls a drone to fly across a target area in a gridded manner, following a preset flight path and altitude. During the drone's flight, the system simultaneously triggers an optical camera to acquire ground optical images at fixed time intervals, forming a raw remote sensing image sequence. Simultaneously, the system triggers a lidar to acquire three-dimensional spatial coordinate points at a fixed pulse frequency, forming a raw laser point cloud data set. The system records the instantaneous position coordinates and attitude angle data of the drone during each acquisition of optical images and laser point clouds. The system then stores the raw remote sensing image sequence and the raw laser point cloud data set in time-aligned format according to the acquisition timestamp. The process of performing stereo pair matching on the original remote sensing image sequence and generating an initial sparse point cloud 3D coordinate set is as follows: Stereo pair sets are formed by extracting overlapping image pairs between adjacent flight strips from the original remote sensing image sequence. Feature point detection is performed on each pair of overlapping images in the stereo pair set to extract scale-invariant feature points. Nearest neighbor matching is performed on the scale-invariant feature points detected between two adjacent images to establish corresponding relationships between corresponding points. The disparity value of each matched point pair is calculated based on the pixel coordinate differences corresponding to the corresponding points. The ground 3D coordinates corresponding to each matched point pair are calculated using the forward intersection method, combining the instantaneous position coordinates and attitude angle data of the UAV with the disparity value. All calculated ground 3D coordinates are then aggregated into an initial sparse point cloud 3D coordinate set.

[0094] In one example scenario, the target survey area is a sparsely vegetated hilly area covering approximately three square kilometers. A drone equipped with a 40-megapixel optical camera and a lidar with a ranging accuracy of two centimeters is used to perform a grid-based, full-coverage flight over the target survey area, following a pre-planned serpentine flight path and a relative flight altitude of 200 meters. In practice, the drone is controlled to perform this grid-based, full-coverage flight over the target survey area according to the preset flight path and altitude. During the drone's flight, the optical camera is simultaneously triggered to acquire ground optical images at fixed time intervals, forming a raw remote sensing image sequence. Simultaneously, the lidar is triggered to acquire three-dimensional spatial coordinate points at a fixed pulse frequency, forming a raw laser point cloud data set. The instantaneous position coordinates and attitude angle data of the drone are recorded each time optical images are acquired; the instantaneous position coordinates and attitude angle data of the drone are also recorded each time laser point cloud data is acquired. The raw remote sensing image sequence and the raw laser point cloud data set are stored with time alignment according to the acquisition timestamp.

[0095] Stereo image pairs are formed by extracting overlapping image pairs between adjacent flight strips from the original remote sensing image sequence. Scale-invariant feature point detection is performed on each pair of overlapping images in the stereo image pair set to extract scale-invariant feature points. Nearest neighbor matching is performed on the scale-invariant feature points detected between two adjacent images to establish corresponding point relationships. The disparity value of each matched point pair is calculated based on the pixel coordinate difference corresponding to the corresponding points. The ground 3D coordinates corresponding to each matched point pair are calculated using forward intersection method by combining the instantaneous position coordinates and attitude angle data of the UAV with the disparity value. In a specific calculation process, the photogrammetric depth value is obtained based on the geometric relationship between disparity and baseline, expressed by the following expression:

[0096]

[0097] in: For photographic depth value, The baseline length components of the left and right projection centers of the stereo image pair in the horizontal flight direction are given. The focal length of an optical camera. To match the disparity values ​​of point pairs, the planar coordinates of ground points in the image space auxiliary coordinate system are determined by the normalized relationship between the photographic depth value and the image point coordinates after rotation by the attitude angle. Then, the image space auxiliary coordinates are transformed to the local projection coordinate system using the rotation matrix composed of the projection center position and attitude angle recorded by the UAV, thus obtaining the ground three-dimensional coordinates corresponding to each matching point pair. All the calculated ground three-dimensional coordinates are then collected into an initial sparse point cloud three-dimensional coordinate set.

[0098] In one embodiment of the present invention, the process of performing a spatial coordinate system-1 transformation between the initial sparse point cloud 3D coordinate set and the original laser point cloud data set is as follows: extracting the geodetic latitude, longitude, and elevation coordinates of each point in the initial sparse point cloud 3D coordinate set; extracting the local projected coordinate system plane coordinates and elevation coordinates of each point in the original laser point cloud data set; selecting at least three common control points between the initial sparse point cloud 3D coordinate set and the original laser point cloud data set; calculating a seven-parameter transformation matrix between the geodetic coordinate system and the local projected coordinate system based on the common control points; transforming all points in the initial sparse point cloud 3D coordinate set to the local projected coordinate system according to the seven-parameter transformation matrix; and superimposing and storing the initial sparse point cloud 3D coordinate set after coordinate system-1 transformation with the original laser point cloud data set.

[0099] The process of extracting surface reflection intensity features and texture gradient features from the data after coordinate system transformation and constructing a joint feature matrix is ​​as follows: Extract the echo reflection intensity value of each laser point from the original laser point cloud data set after coordinate system transformation, and map the echo reflection intensity value of each laser point onto a two-dimensional regular grid according to the planar coordinate position of the point to form a reflection intensity feature map; extract the original remote sensing image pixel value corresponding to each point from the initial sparse point cloud three-dimensional coordinate set after coordinate system transformation, and calculate the spatial change rate of each pixel value in the horizontal and vertical directions as the texture gradient feature; concatenate the reflection intensity feature value and the texture gradient feature value at the same grid position according to the grid row and column order to form a two-dimensional feature vector group, and organize the two-dimensional feature vector groups of all grids into a three-dimensional tensor according to the spatial arrangement order as the joint feature matrix.

[0100] Continuing with the example scenario of the hilly survey area, after the UAV completes data collection, it obtains an initial sparse point cloud 3D coordinate set with geodetic latitude, longitude, and elevation coordinates, and an original laser point cloud data set with local projected coordinate system plane coordinates and elevation values. Since the two sets of point cloud data are located in different spatial reference frames, a spatial coordinate system transformation is performed to align the geometric reference of the data. In specific implementation, the geodetic latitude, longitude, and elevation coordinates of each point in the initial sparse point cloud 3D coordinate set are extracted, and the local projected coordinate system plane coordinates and elevation coordinates of each point in the original laser point cloud data set are extracted. At least three common control points between the initial sparse point cloud 3D coordinate set and the original laser point cloud data set are selected. These common control points are selected from ground markers that have been pre-deployed and accurately measured in both coordinate systems within the survey area. A seven-parameter transformation matrix between the geodetic coordinate system and the local projected coordinate system is calculated based on the common control points. All points in the initial sparse point cloud 3D coordinate set are transformed to the local projected coordinate system according to the seven-parameter transformation matrix. The initial sparse point cloud 3D coordinate set after coordinate system transformation is superimposed and stored with the original laser point cloud data set.

[0101] The echo reflection intensity value of each laser point is extracted from the original laser point cloud data set after coordinate system transformation. This value is then mapped onto a two-dimensional regular grid according to the corresponding planar coordinates of the laser point to form a reflection intensity feature map. In a specific feature map generation process, the resolution of the two-dimensional regular grid is set to match the preset grid spacing of the subsequent raster-based digital elevation model. When multiple laser points exist within a grid location, the arithmetic mean of the echo reflection intensity values ​​of these multiple laser points is taken as the reflection intensity feature value for that grid location. When no laser point exists within a grid location, the reflection intensity feature value of the blank grid is filled using inverse distance weighted interpolation of neighboring grids with values. The original remote sensing image pixel value corresponding to each point is extracted from the initial sparse point cloud three-dimensional coordinate set after coordinate system transformation. The spatial variation rate of each pixel value in the horizontal and vertical directions is calculated as a texture gradient feature. The calculation of the horizontal and vertical spatial variation rates is expressed by the following expressions:

[0102]

[0103] in: The magnitude of the rate of change in the horizontal direction. The magnitude of the rate of change in the vertical direction. and The grid positions are respectively The pixel values ​​of adjacent pixels on the right and left sides. and The grid positions are respectively The pixel values ​​of the adjacent pixels below and above will and The two-dimensional vector synthesized into texture gradient features is used as the texture gradient feature value of the corresponding grid position; the reflection intensity feature value and the texture gradient feature value at the same grid position are concatenated according to the grid row and column order to form a two-dimensional feature vector group; and the two-dimensional feature vector groups of all grids are organized into a three-dimensional tensor according to the spatial arrangement order as the joint feature matrix.

[0104] In one embodiment of the present invention, the process of inputting the joint feature matrix into a pre-trained deep convolutional neural network and outputting a preliminary terrain category probability distribution map is as follows: A deep convolutional neural network model pre-trained on a labeled terrain sample dataset is obtained. The training process of this model includes feature extraction and enhancement of the joint feature matrix through a multi-scale feature fusion module. The joint feature matrix is ​​cut into multiple fixed-size feature patches according to the spatial coverage of the original remote sensing image sequence. Each feature patch is sequentially input into the deep convolutional neural network model for forward propagation calculation. The deep convolutional neural network model outputs a probability vector for each feature patch, where each component corresponds to a preset terrain category. The output probability vectors of all feature patches are reassembled into a multi-channel probability map according to their original spatial locations. The terrain category with the highest probability at each spatial location in the multi-channel probability map is marked as the preliminary terrain category at that location, forming a preliminary terrain category probability distribution map.

[0105] The process of separating ground point clouds from non-ground point clouds by performing regional filtering on the original laser point cloud data based on the preliminary terrain category probability distribution map is as follows: (See...) Figure 3 Each laser point in the original laser point cloud data set is located to its corresponding position on the preliminary terrain category probability distribution map according to its planar coordinates. The terrain category with the highest probability at that position is read as the terrain category label of the laser point. When the terrain category label is bare ground or low vegetation, the laser point is marked as a candidate ground point. When the terrain category label is a building or tall vegetation, the laser point is marked as a candidate non-ground point. Progressive morphological filtering is performed on the laser point cloud marked as candidate ground points to remove outliers. Clustering segmentation filtering is performed on the laser point cloud marked as candidate non-ground points to separate independent non-ground object point clouds. The remaining ground points and non-ground points after filtering are stored in the ground point cloud set and the non-ground point cloud set, respectively.

[0106] In the example scenario of the hilly area described above, the joint feature matrix generated in the preceding steps is used to drive terrain classification and subsequent point cloud filtering. This joint feature matrix is ​​a three-dimensional tensor formed by concatenating reflection intensity features and texture gradient features, covering the entire target area. In specific implementation, a deep convolutional neural network model pre-trained on a labeled terrain sample dataset is obtained. The training process of this deep convolutional neural network model includes feature extraction and enhancement of the joint feature matrix through a multi-scale feature fusion module. This multi-scale feature fusion module consists of parallel convolutional layers with different dilation rates, capable of simultaneously capturing both small-scale terrain details and large-scale terrain context features. The joint feature matrix is ​​then divided into multiple fixed-size feature patches according to the spatial coverage of the original remote sensing image sequence. Adjacent feature patches maintain an overlap of one-quarter of the patch's side length at the boundaries to avoid edge effects. Each feature patch is then sequentially input into the deep convolutional neural network model. In the model, forward propagation calculations are performed. The deep convolutional neural network model outputs a probability vector for each feature patch. Each component of the probability vector corresponds to the probability value of the feature patch belonging to a preset terrain category. The preset terrain categories include at least bare ground, low vegetation, buildings, and tall vegetation. The output probability vectors of all feature patches are reassembled into a multi-channel probability map according to their original spatial locations. The probability value of the overlapping area is taken as the arithmetic mean of the output probabilities of multiple feature patches. The terrain category with the highest probability at each spatial location in the multi-channel probability map is marked as the preliminary terrain category for that spatial location. The preliminary terrain category markings for all spatial locations form a preliminary terrain category probability distribution map.

[0107] Based on the preliminary terrain category probability distribution map, the original laser point cloud data set is subjected to regional filtering to separate ground point clouds from non-ground point clouds. In specific implementation, each laser point in the original laser point cloud data set is located to its corresponding spatial position on the preliminary terrain category probability distribution map according to its planar coordinates. The terrain category with the highest probability at the corresponding spatial position is read as the terrain category label of the corresponding laser point. When the terrain category label of a laser point is bare surface or low vegetation, the corresponding laser point is marked as a candidate ground point. When the terrain category label of a laser point is building or tall vegetation, the corresponding laser point is marked as a candidate non-ground point. Progressive morphological filtering is performed on the laser point cloud marked as candidate ground points to eliminate them. Except for outliers, progressive morphological filtering repeatedly filters out outliers above the local terrain trend surface by gradually increasing the filter window size and residual threshold. Clustering segmentation filtering is performed on the laser point cloud marked as candidate non-ground points to separate independent non-ground object point clouds. The clustering segmentation filtering establishes connectivity based on the three-dimensional Euclidean distance between points, dividing points that are more than a preset distance threshold into different clusters. The real ground points retained after progressive morphological filtering are stored in the ground point cloud set, and each cluster obtained by clustering segmentation filtering is stored as an independent non-ground object point cloud in the non-ground point cloud set. The dividing line between the two types of point cloud sets is strictly aligned with the boundary between the exposed surface category and the building category in the preliminary terrain category probability distribution map.

[0108] In one embodiment of the present invention, the process of performing irregular triangular mesh encryption interpolation on the separated ground point cloud and generating a first version of the raster-format digital elevation model is as follows: extracting the planar coordinates and elevation values ​​of all ground points from the ground point cloud set, and constructing a Delaunay irregular triangular mesh using the planar coordinates of all ground points as triangle vertices; performing linear interpolation within each triangle of the Delaunay irregular triangular mesh according to the elevation values ​​of the three vertices of the triangle; arranging regular grid points within the target survey area according to a preset grid spacing, locating the irregular triangular mesh triangle in which each regular grid point is located, and calculating the interpolated elevation value of the triangle; organizing the elevation values ​​of all regular grid points into a two-dimensional matrix according to the grid row and column order; performing inverse distance weighted interpolation on the blank grid points in the two-dimensional matrix to fill the blanks, and using the two-dimensional matrix after filling the blanks as the first version of the raster-format digital elevation model.

[0109] The process of calculating the elevation residual distribution field between the first version of the raster digital elevation model and the initial sparse point cloud 3D coordinate set is as follows: The elevation value of each regular grid point in the first version of the raster digital elevation model is used as the reference elevation value, and the elevation value of each sparse point in the initial sparse point cloud 3D coordinate set is used as the observed elevation value. For each sparse point, the grid point to which its planar coordinates belong in the first version of the raster digital elevation model is located is found, and the difference between the observed elevation value and the reference elevation value of the sparse point is calculated as the elevation residual of that point. The elevation residuals of all sparse points are assigned to corresponding grid cells according to their planar coordinate positions. The arithmetic mean of multiple elevation residuals assigned to the same grid cell is calculated as the representative residual value of that grid cell. Grid cells with representative residual values ​​are combined with null grid cells without sparse points to form an incomplete residual field. Kriging space interpolation is performed on the null grid cells in the incomplete residual field to generate a complete gridded residual field as the elevation residual distribution field.

[0110] In the example scenario of the hilly area described above, the ground point cloud set has been separated by regional filtering and used to construct the first version of the raster-based digital elevation model. The initial sparse point cloud 3D coordinate set also participates in the calculation of the elevation residual distribution field. In specific implementation, the planar coordinates and elevation values ​​of all ground points are extracted from the ground point cloud set, and the planar coordinates of all ground points are used as the vertices of triangles to construct a Delaunay irregular triangular mesh. The circumcircle of each triangle in the Delaunay irregular triangular mesh does not contain other ground points, ensuring the best fit of the triangular mesh to the planar distribution of ground points. Linear interpolation is performed inside each triangle of the Delaunay irregular triangular mesh according to the elevation values ​​of the three vertices of the corresponding triangle. The linear interpolation uses the inclined plane formed by the elevations of the three vertices of the triangle to determine the elevation of any position inside the triangle. Within the target survey area, the grid is arranged according to the preset grid spacing rules. The grid points are set with a preset grid spacing of one meter. The planar coordinates of each regular grid point are determined by accumulating the horizontal and vertical coordinates along the grid spacing from the origin at the lower left corner of the target area. For each regular grid point, the Delaunay irregular triangulation triangle it belongs to is located. The grid point is assigned to a triangle by determining whether its planar coordinates fall inside the convex hull formed by the three vertices of the triangle. After location, the elevation value of the regular grid point is calculated by using the linear interpolation relationship of the corresponding triangle. The elevation values ​​of all regular grid points are organized into a two-dimensional matrix according to the grid row and column order. The row index of the two-dimensional matrix corresponds to the north-south direction, and the column index corresponds to the east-west direction.

[0111] Optionally, after organizing the elevation values ​​of all regular grid points into a two-dimensional matrix, blank grid points in the two-dimensional matrix that are not covered by the Delaunay irregular triangular mesh at the edge corners of the target survey area need to be filled by inverse distance weighted interpolation. For each blank grid point, inverse distance weighted interpolation searches for existing elevation value grid points within a certain radius around it as interpolation reference points, and calculates the elevation value of the blank grid point using a power-weighted distance method. The calculation of inverse distance weighted interpolation is expressed by the following expression:

[0112]

[0113] in: blank grid dots interpolated elevation values, This represents the total number of existing elevation value grid points found within the search radius. For the first The elevation values ​​of existing elevation grid points. blank grid dots With the The planar Euclidean distance between existing elevation grid points. The power constant for the distance weight is set to two; the two-dimensional matrix after filling the blank grid points is used as the first version of the raster digital elevation model.

[0114] In practical implementation, the elevation value of each regular grid point in the first version of the raster-style digital elevation model is used as the reference elevation value, and the elevation value of each sparse point in the initial sparse point cloud 3D coordinate set is used as the observed elevation value. For each sparse point in the initial sparse point cloud 3D coordinate set, the sparse point's planar coordinates are used to locate it in the grid array of the first version of the raster-style digital elevation model. The coverage area of ​​which regular grid point the sparse point's planar coordinates fall within is determined, and the corresponding regular grid point is identified as the grid point to which the sparse point belongs. The difference between the observed elevation value of the sparse point and the reference elevation value of its grid point is calculated, and this difference is used as the elevation residual of the corresponding sparse point. The elevation residuals of all sparse points are distributed to the corresponding grid cells according to the planar coordinate positions of the sparse points. The division of the grid cells is completely consistent with the grid division of the first version of the raster-style digital elevation model. Each grid cell... A cell may contain multiple sparse points or no sparse points. The arithmetic mean of multiple elevation residuals assigned to the same grid cell is calculated, and this arithmetic mean is used as the representative residual value for the corresponding grid cell. Grid cells with representative residual values ​​are combined with null-value grid cells without sparse points to form an incomplete residual field. The incomplete residual field is a two-dimensional residual matrix with the same number of rows and columns as the first version of the raster-based digital elevation model. Grid cells with values ​​represent the representative residual values ​​of their respective grids, while grid cells without values ​​represent null values. Kriging space interpolation is performed on the null-value grid cells in the incomplete residual field. Based on the known spatial distribution structure of the representative residual values ​​and the semi-variogram model, Kriging space interpolation performs optimal linear unbiased estimation for each null-value grid cell. After interpolation, all null-value grid cells are assigned estimated residual values, thus generating a complete gridded residual field as the elevation residual distribution field.

[0115] In one embodiment of the present invention, the process of using the elevation residual distribution field to correct the elevation values ​​of the first raster digital elevation model grid by grid to generate the second raster digital elevation model, and outputting the second raster digital elevation model as the final terrain elevation data inversion calculation result, is as follows: The residual value of each grid cell in the elevation residual distribution field is read from that grid cell; the reference elevation value of the same grid cell in the first raster digital elevation model is algebraically summed with the residual value of that grid cell; the summation result is written as the corrected elevation value of that grid cell into the corresponding position of a new blank two-dimensional matrix; the steps of reading the residual value and algebraically summing are repeated for each grid cell in the blank two-dimensional matrix until all grid cells have been processed; the two-dimensional matrix containing all the corrected elevation values ​​is used as the second raster digital elevation model; the second raster digital elevation model is exported as a file according to the standard geographic information system raster format; and the exported file is output as the final terrain elevation data inversion calculation result.

[0116] In the example scenario of the hilly area mentioned above, the first version of the raster digital elevation model and the elevation residual distribution field have been constructed. Both are two-dimensional matrices with the same number of rows and columns. Each residual value in the elevation residual distribution field quantitatively represents the systematic elevation deviation of the first version of the raster digital elevation model relative to the initial sparse point cloud three-dimensional coordinate set at the corresponding spatial location. The elevation residual distribution field is used to correct the elevation values ​​of the first version of the raster digital elevation model grid by grid.

[0117] In practical implementation, the residual value of each grid cell in the elevation residual distribution field is read from the two-dimensional matrix of the elevation residual distribution field according to the row and column indices. The reference elevation value of the grid cell with the same row and column index in the first version of the raster digital elevation model is algebraically summed with the residual value of the corresponding grid cell. The algebraic summation is expressed by the following expression:

[0118]

[0119] in: In the revised second edition of the raster digital elevation model, the first... Line 1 Corrected elevation values ​​for grid cells. The first version of the raster digital elevation model Line 1 Reference elevation values ​​for grid cells, The first in the elevation residual distribution field Line 1 The residual value of the grid cell. and These are the row and column indices, respectively. The corrected elevation values ​​obtained by algebraic summation are written into the corresponding positions of a new blank two-dimensional matrix. The row and column dimensions of the blank two-dimensional matrix are exactly the same as those of the first version of the raster digital elevation model. For each grid cell in the blank two-dimensional matrix, the steps of reading the residual value, performing algebraic summation, and writing the corrected elevation value are repeated until all grid cells are processed. During the processing, the process is traversed row by row, first fixing the row index and incrementing the column index to the end of the row before jumping to the next row. The two-dimensional matrix with all the corrected elevation values ​​written into it is used as the second version of the raster digital elevation model.

[0120] Optionally, after completing the grid-by-grid correction, the second version of the raster digital elevation model is subjected to edge smoothing. When the elevation difference between adjacent grid cells exceeds the preset allowable threshold for terrain change, the elevation value of the central grid cell is replaced with the moving weighted average of the elevation values ​​of the adjacent grid cells to eliminate local micro-undulations that may be caused by residual field interpolation.

[0121] The second version of the raster digital elevation model is exported as a file according to the standard geographic information system raster format. The header information of the exported raster file includes the number of raster rows and columns, the side length of the raster cell, the coordinates of the starting plane at the lower left corner, the identifier of no data value, and the description of the elevation datum. The raster data is written into the corrected elevation value of each grid cell in row order. The exported file is output as the final topographic elevation data inversion calculation result.

[0122] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments that can be applied to other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for inverting and calculating terrain elevation data assisted by UAV remote sensing images, characterized in that, include: The original remote sensing image sequence and original laser point cloud data set captured by the optical camera and lidar on the UAV in the target survey area are collected. The original remote sensing image sequence is subjected to stereo image pair matching and an initial sparse point cloud three-dimensional coordinate set is generated. The initial sparse point cloud 3D coordinate set is transformed with the original laser point cloud data set into a spatial coordinate system. The surface reflection intensity features and texture gradient features of the ground objects are extracted from the data after the coordinate system transformation, and a joint feature matrix is ​​constructed. The joint feature matrix is ​​input into a pre-trained deep convolutional neural network and a preliminary terrain category probability distribution map is output. Based on the preliminary terrain category probability distribution map, the original laser point cloud data set is filtered by region to separate ground point clouds from non-ground point clouds; The separated ground point cloud was subjected to irregular triangular mesh encryption interpolation to generate the first version of the raster digital elevation model; Calculate the elevation residual distribution field between the first version of the raster-format digital elevation model and the initial sparse point cloud three-dimensional coordinate set; The first raster digital elevation model is corrected grid by grid using the elevation residual distribution field to generate a second raster digital elevation model. The second raster digital elevation model is then output as the final terrain elevation data inversion calculation result.

2. The method for inverting terrain elevation data assisted by UAV remote sensing images according to claim 1, characterized in that, The steps for collecting raw remote sensing image sequences and raw laser point cloud data sets captured by the optical camera and lidar onboard a drone on the target survey area specifically include: Control the drone to perform grid-based full-coverage flight over the target survey area according to the preset route and preset flight altitude; During the flight of the UAV, the optical camera is synchronously triggered to acquire ground optical images at fixed time intervals and form the original remote sensing image sequence. During the flight of the UAV, the lidar is synchronously triggered to collect three-dimensional spatial coordinate points at a fixed pulse frequency and form the original laser point cloud data set. Record the instantaneous position coordinates and attitude angle data of the UAV each time optical images are acquired; Record the instantaneous position coordinates and attitude angle data of the UAV each time laser point cloud is acquired; The original remote sensing image sequence and the original laser point cloud data set are stored in time-aligned according to the acquisition timestamp.

3. The method for inverting terrain elevation data assisted by UAV remote sensing images according to claim 2, characterized in that, The steps of performing stereo pair matching on the original remote sensing image sequence and generating an initial sparse point cloud 3D coordinate set specifically include: Extract overlapping image pairs between adjacent flight strips from the original remote sensing image sequence to form a stereo image pair set; Feature point detection is performed on each pair of overlapping images in the stereo image pair set to extract scale-invariant feature points; Nearest neighbor matching is performed on scale-invariant feature points detected between two adjacent images to establish correspondences of corresponding points. The disparity value of each matching point pair is calculated based on the difference in pixel coordinates corresponding to the same point. The ground 3D coordinates corresponding to each matching point are calculated by combining the instantaneous position coordinates and attitude angle data of the UAV with the parallax value using the forward intersection method; All the calculated 3D ground coordinates are aggregated into the initial sparse point cloud 3D coordinate set.

4. The method for inverting terrain elevation data assisted by UAV remote sensing images according to claim 1, characterized in that, The step of performing a spatial coordinate system transformation between the initial sparse point cloud 3D coordinate set and the original laser point cloud data set specifically includes: Extract the latitude, longitude, and elevation coordinates of each point in the initial sparse point cloud 3D coordinate set; Extract the local projected coordinate system plane coordinates and elevation coordinates of each point in the original laser point cloud data set; Select at least three common control points between the initial sparse point cloud 3D coordinate set and the original laser point cloud data set; Calculate the seven-parameter transformation matrix between the geodetic coordinate system and the local projected coordinate system based on the aforementioned common control points; Transform all points in the initial sparse point cloud 3D coordinate set to the local projected coordinate system according to the seven-parameter transformation matrix; The initial sparse point cloud 3D coordinate set after coordinate system transformation is superimposed and stored with the original laser point cloud data set.

5. The method for inverting terrain elevation data assisted by UAV remote sensing images according to claim 4, characterized in that, The steps for extracting surface reflectance intensity features and texture gradient features from data after coordinate system transformation and constructing a joint feature matrix specifically include: Extract the echo reflection intensity value of each laser point from the original laser point cloud data set after coordinate system transformation; The echo reflection intensity value of each laser point is mapped onto a two-dimensional regular grid according to the planar coordinate position of that point to form a reflection intensity feature map; Extract the original remote sensing image pixel value corresponding to each point from the initial sparse point cloud three-dimensional coordinate set after coordinate system transformation; Calculate the spatial rate of change of each pixel value in the horizontal and vertical directions as texture gradient features; The reflection intensity feature value and the texture gradient feature value at the same grid location are concatenated in the order of grid rows and columns to form a two-dimensional feature vector group; The two-dimensional feature vectors of all grids are arranged in spatial order to form a three-dimensional tensor as the joint feature matrix.

6. The method for inverting terrain elevation data assisted by UAV remote sensing images according to claim 5, characterized in that, The steps of inputting the joint feature matrix into a pre-trained deep convolutional neural network and outputting a preliminary terrain category probability distribution map specifically include: Obtain a deep convolutional neural network model that has been pre-trained on a labeled terrain sample dataset; The joint feature matrix is ​​cut into multiple fixed-size feature patches according to the spatial coverage of the original remote sensing image sequence; Each feature map block is sequentially input into the deep convolutional neural network model for forward propagation calculation; The deep convolutional neural network model outputs a probability vector for each feature map block; Each component in the probability vector corresponds to a feature patch belonging to a preset terrain category; The output probability vectors of all feature map blocks are reassembled into a multi-channel probability map according to their original spatial positions. The terrain category with the highest probability at each spatial location in the multi-channel probability map is marked as the preliminary terrain category at that location, and the preliminary terrain category probability distribution map is formed. The training process of the deep convolutional neural network model includes: extracting and enhancing features from the joint feature matrix through a multi-scale feature fusion module.

7. The method for inverting terrain elevation data assisted by UAV remote sensing images according to claim 6, characterized in that, The step of performing regional filtering on the original laser point cloud data group based on the preliminary terrain category probability distribution map to separate ground point clouds from non-ground point clouds specifically includes: Each laser point in the original laser point cloud data set is located at its corresponding position in the preliminary terrain category probability distribution map according to its planar coordinates; Read the terrain category with the highest probability at that location and use it as the terrain category label for that laser point; When the terrain category label is bare ground or low vegetation, the laser point is marked as a candidate ground point; When the terrain category label is building or tall vegetation, the laser point is marked as a candidate non-ground point; Progressive morphological filtering is performed on the laser point cloud marked as candidate ground points to remove outliers; Clustering segmentation filtering is performed on the laser point cloud labeled as candidate non-ground points to separate the independent non-ground object point cloud; The remaining ground points and non-ground points after filtering are stored in the ground point cloud set and the non-ground point cloud set, respectively.

8. The method for inverting terrain elevation data assisted by UAV remote sensing images according to claim 7, characterized in that, The steps for performing irregular triangular mesh encryption and interpolation on the separated ground point cloud to generate the first version of the raster digital elevation model specifically include: Extract the planar coordinates and elevation values ​​of all ground points from the ground point cloud set; Construct a Delaunay irregular triangulation using the planar coordinates of all ground points as the vertices of triangles; Linear interpolation is performed inside each triangle of the Delaunay irregular triangular network according to the elevation values ​​of the three vertices of the triangle; Arrange regular grid points within the target survey area according to the preset grid spacing; For each regular grid point, locate the irregular triangular network triangle it belongs to and calculate the interpolated elevation value inside the triangle; Organize the elevation values ​​of all regular grid points into a two-dimensional matrix according to the grid row and column order; Inverse distance weighted interpolation is performed on the blank grid points in the two-dimensional matrix to fill the blanks. The two-dimensional matrix after filling the blanks is used as the first version of the raster digital elevation model.

9. The method for inverting terrain elevation data assisted by UAV remote sensing images according to claim 8, characterized in that, The steps for calculating the elevation residual distribution field between the first raster-format digital elevation model and the initial sparse point cloud 3D coordinate set specifically include: The elevation value of each regular grid point in the first version of the raster digital elevation model is used as the reference elevation value. The elevation value of each sparse point in the initial sparse point cloud three-dimensional coordinate set is used as the observed elevation value. For each sparse point, find the grid point in the first version of the raster digital elevation model to which the plane coordinates of that sparse point belong; The difference between the observed elevation value of the sparse point and the reference elevation value of the grid point is calculated as the elevation residual of the point. Allocate the elevation residuals of all sparse points to the corresponding grid cells according to their planar coordinate positions; The arithmetic mean of multiple elevation residuals assigned to the same grid cell is calculated and used as the representative residual value of that grid cell; An incomplete residual field is formed by combining grid cells that represent residual values ​​with null grid cells that do not have sparse points. Kriging spatial interpolation is performed on the null grid cells in the incomplete residual field to generate a complete gridded residual field as the elevation residual distribution field.

10. The method for inverting terrain elevation data assisted by UAV remote sensing images according to claim 9, characterized in that, The steps of using the elevation residual distribution field to perform grid-by-grid elevation value correction on the first raster digital elevation model to generate a second raster digital elevation model, and outputting the second raster digital elevation model as the final terrain elevation data inversion calculation result, specifically include: The residual value of each grid cell in the elevation residual distribution field is read from that grid cell; The reference elevation value of the same grid cell in the first version of the raster digital elevation model is algebraically summed with the residual value of that grid cell. The summation result is written as the corrected elevation value of the grid cell into the corresponding position of a new blank two-dimensional matrix; The steps of reading the residual value and calculating the algebraic summation are repeated for each grid cell in the blank two-dimensional matrix until all grid cells have been processed. The two-dimensional matrix containing all the corrected elevation values ​​is used as the second version of the raster-format digital elevation model. Export the second version of the raster digital elevation model as a file according to the standard geographic information system raster format; The exported file is output as the final terrain elevation data inversion calculation result.

Images