Hyperspectral image-point cloud stereo registration method based on light ray tracing correction

Abstract

The application discloses a hyperspectral image-point cloud stereo registration method based on light ray tracking correction, and belongs to the technical field of stereo registration of hyperspectral images and laser radar data. In order to solve the problem that stereo registration of hyperspectral images and laser radar point clouds of a low-altitude airborne platform has a big loss of stereo dimension information caused by low-altitude visual angle error, the application first calculates the spatial position and flight direction vector of the unmanned aerial vehicle when scanning a single frame of hyperspectral image, determines a layered plane according to a plane normal vector and a fixed point on the plane, and determines a spatial layer corresponding to a single hyperspectral scanning frame based on the layered plane. Then, according to any point in the hyperspectral image point cloud layer, a corresponding offset point is obtained on the layered result, and a corresponding relationship between each point and a grid is established by orthographic projection according to coordinate grid parameters, and then a digital surface model is obtained. According to the phase consistency principle, edge extraction is performed on the digital surface model and the hyperspectral image, and then registration is performed on the digital surface model and the hyperspectral image.

Description

Technical Field

[0001] This invention belongs to the field of stereo registration technology of hyperspectral images and lidar data, specifically relating to a stereo registration method for hyperspectral images and point clouds based on ray tracing correction. Background Technology

[0002] Hyperspectral imaging can acquire two-dimensional and spectral information of the observed scene / target, resulting in a unified remote sensing image. Lidar, a technology that uses a laser beam to scan a target object and receives the reflected laser signal to quickly acquire three-dimensional spatial information of the observed scene / target, offers advantages such as high precision and immunity to weather conditions. To leverage the respective strengths of hyperspectral imaging and lidar, the joint processing of hyperspectral images and lidar data is of significant research importance, capable of generating both spectral and three-dimensional spatial information of the observed scene / target, and holds broad application prospects.

[0003] Because hyperspectral sensors and lidar operate on different principles, registration is required before joint information processing. Current hyperspectral image-lidar data registration primarily focuses on two-dimensional registration for spaceborne / high-altitude airborne platforms, with relatively little research on stereo registration for low-altitude airborne platforms. Low-altitude airborne hyperspectral-lidar can acquire higher spatial resolution spatial three-dimensional and spectral information, which is significant in fields such as 3D reconstruction and digital twins. However, existing research on stereo registration of hyperspectral-lidar data for low-altitude UAV platforms is still insufficient. Directly using mainstream image registration methods struggles to overcome low-altitude viewpoint errors, leading to loss of stereo dimensional information. Furthermore, mainstream stereo registration methods still rely heavily on manual calibration-mapping processes, requiring manual intervention and unsuitable for real-time data registration. Summary of the Invention

[0004] This invention addresses the problem of significant loss of stereoscopic dimensional information due to low-altitude viewing angle errors during stereo registration of hyperspectral images and lidar point clouds from low-altitude airborne platforms.

[0005] A hyperspectral image-point cloud stereo registration method based on ray tracing correction includes the following steps:

[0006] S1. Layer the lidar point cloud, including the following steps:

[0007] S1.1 Calculate the spatial position and flight direction vector of the UAV when scanning a single hyperspectral frame:

[0008] Let the trajectory data be {P} m ,T m POS m}, where m = 1, 2, ..., M, M represents the number of hyperspectral image frames; P m =(Em N m H m () represents the coordinates of the E, N, H axes in the UTM coordinate system; T m For timestamps; ω is the direction angle of flight. m ,κ m , These are pitch angle, heading angle, and yaw angle, respectively.

[0009] The direction vector at each frame scan time is calculated using the following formula:

[0010]

[0011] Among them, v E ,v N ,v H Represents the unit direction vectors of the E, N, and H axes in UTM coordinates;

[0012] S1.2, Spatial layer Z corresponding to a single hyperspectral scan frame m Represented as front and rear layered planes Π m Π m+1 Intermediate mezzanine space, based on layered plane Π i Determine the spatial layer Z corresponding to a single hyperspectral scan frame m ;

[0013] Layered plane Π i From the plane normal vector d i With a fixed point P on the plane i Determine, where i = 1, 2, ..., M, M+1;

[0014] S2. Based on the layered results, generate a digital surface model based on the offset projection of ray-traced point clouds and establish a correspondence with the original point cloud. The process of generating a digital surface model based on the offset projection of ray-traced point clouds includes the following steps:

[0015] S2.1. Based on the point cloud layering results from step 1, and based on any point p within the point cloud layer corresponding to a certain frame of the hyperspectral image... n =(E n N n H n ) Obtain the corresponding offset point Calculate using the following formula:

[0016]

[0017] Where κ' is the track angle, representing the angle relative to true north. N The acute angle formed; H pj P represents the height of the hyperspectral imaging plane. m =(Em N m H m () represents the spatial coordinates of the UAV at the moment of scanning this frame;

[0018] Iterate through all point clouds within the layer, perform offset transformations on these point clouds, and preserve their height H. n Without changing anything, generate an offset point cloud;

[0019] S2.2. Establish a raster with R×C rows and columns, where R×C should be greater than the resolution of the hyperspectral image; based on the coordinate raster parameters, establish the correspondence between each point and the raster using orthophoto projection of the offset point cloud: in, Indicates lidar points The coordinates of the corresponding grid in the digital surface model; (r,c) represents the horizontal and vertical coordinates of a certain grid in the lidar digital surface model;

[0020] The value of a single grid cell is the highest value among all the height values ​​of all points within the cell, which is denoted as the gray level of the grid cell. After traversing all grid cells, the gray level is linearly mapped to [0,1] to obtain the digital surface model.

[0021] S3. Based on the principle of phase consistency, edge extraction is performed on the digital surface model and the hyperspectral image;

[0022] S4. Using image registration methods, the digital surface model is registered with the hyperspectral image to generate a hyperspectral point cloud.

[0023] Furthermore, in S1.2, the layered plane Π i normal vector d on the plane i Determined in the following manner:

[0024] Layered plane Π2 to Π M The corresponding plane normal vector is the flight direction vector d at the previous and next time steps. m The average vector;

[0025] First and last planes Π1 and Π M+1 The normal vector is taken as the flight direction vectors d1 and d2 at the beginning and end of the flight time. M .

[0026] Furthermore, in S1.2, the layered plane Π i Fixed point P on i Determined in the following manner:

[0027] Layered plane Π2 to Π M The fixed point on the plane is the flight coordinate P at the previous and next time points. m The midpoint;

[0028] The fixed point of the first plane Π1 is determined by extending the fixed points of Π2 and Π3 in the opposite direction and taking points of equal length;

[0029] End plane Π M+1 The fixed point of the surface is Π M-1 Π M The fixed point is extended in the opposite direction to find a point of equal length to determine the point.

[0030] Furthermore, the layered plane Π described in S1.2 i From the plane normal vector d i With a fixed point P on the plane i The determination process includes the following steps:

[0031] Based on the determined plane normal vector d i Using the plane normal vector d i =A i ·v E +B i ·v N +C i ·v H Obtain coefficient A i B i C i Based on a fixed point P i =(E i N i H i ), (E i N i H i The coordinates of the E, N, and H axes in the UTM coordinate system, and the plane Π i as follows:

[0032] Π i :π i =A i ·e+B i ·n+C i ·h-(E i +N i +H i ) = 0

[0033] Where, π i =A i ·e+B i ·n+C i ·h-(E i +N i +H i ) represents the plane Π i The corresponding function, π i (p) represents substituting the coordinates of a point p = (E, N, H) in space into the function π. iThe calculated function values; e, n, and h are the coordinates of the E, N, and H axes in the three UTM coordinate system, expressed as independent variables in the function form.

[0034] Based on layered plane Π i Determine the spatial layer Z corresponding to a single hyperspectral scan frame m The process includes the following steps:

[0035] Based on M+1 boundary planes, the point cloud space is divided into M scanning layers Z corresponding to the hyperspectral image frames. m Each scan layer consists of a plane Π m With Π m+1 Determine any point p in the point cloud. n =(E n N n H n In the scanning layer Z m The necessary and sufficient condition within is expressed as:

[0036]

[0037] in, This indicates that the two are equivalent;

[0038] Based on the order of frame scanning, the corresponding formulas for each point in the spatial layer are substituted for verification. For points p that satisfy the formulas... n Record the corresponding scanning layer number, which is not used for subsequent layer verification; calculate the correspondence between each lidar point cloud and the spatial layer based on the layer plane information.

[0039] or,

[0040] In step 1.2, if the UAV flight satisfies the straight-line flight attitude stability conditions (a) and (b), the fast point cloud layering algorithm is used to determine the layering plane Π. i And based on the hierarchical plane Π i Determine the spatial layer Z corresponding to a single hyperspectral scan frame m ;

[0041] The drone's flight meets the following attitude stability conditions for straight-line flight:

[0042] (a) The flight trajectory can approximate a straight path, that is, it satisfies the instantaneous heading angle κ of all trajectory points in the path. m It tends towards a constant K, with a deviation not exceeding ±5°;

[0043] (b) During a single-segment straight flight, the altitude fluctuation does not exceed 5 meters, and the flight pitch angle ω is constant at all scanning times. m Approaching 0, and all not exceeding 10°;

[0044] Determine the layering plane Π using a fast point cloud layering algorithmi The process includes the following steps:

[0045] When the flight path is relatively stable and approximately straight, at this point, the layered planes are perpendicular to the ground surface, and all planes are parallel to each other, i.e., C i =0, A i B i Let A and B be constants, then plane Π i From fixed point P i Confirmed, indicated as:

[0046] Π i :π i =A·e+B·nD i =0

[0047] Among them, the constant term D of the plane equation i =E i +N i +H i Let E, N, and H be the coordinates of fixed points, and let A = sinK and B = cosK.

[0048] Based on layered plane Π i Determine the spatial layer Z corresponding to a single hyperspectral scan frame m Includes the following steps:

[0049] During straight-line flight, the plane group constant term D i It is a monotonic sequence; it passes through any point p in space. n And the plane parallel to the plane group is represented as vtp n :A·e+B·nD n =0; where the constant term D n Let E, N, and H be the coordinates of this point.

[0050] At this time p n In scan layer Z m The necessary and sufficient condition within is expressed as:

[0051]

[0052] Among them, D m D m+1 For scanning layer Z m The corresponding front and rear layer plane equations Π m Π m+1 constant term;

[0053] When the hyperspectral imager flight data meets the straight-line attitude stability condition, the point cloud is first traversed to obtain preliminary point cloud layering results p. n ∈Z m And then according to The non-approximate plane equations in the model are used to verify the correctness of the layering results; if the verification is incorrect, then p... n Substitute the nearest neighbor layer Z m±t Continue until it is verified to be correct.

[0054] Furthermore, the process of edge extraction from the digital surface model and hyperspectral image includes the following steps:

[0055] Single-band images in hyperspectral images With digital surface model I DSM All images are grayscale images, where S represents the total number of bands in the hyperspectral image. A two-dimensional Log-Gabor filter is used to filter these grayscale images, and an edge map M is generated using a phase consistency algorithm. The edge map group of the hyperspectral image is denoted as... The edge graph of the digital surface model is M DSM ;

[0056] Let the noise threshold group of the highlight image edge group be . Noise suppression is achieved by iterating through the entire band edge map of the highlight image using the following formula:

[0057]

[0058] in, This is the edge map of the hyperspectral image after full-band noise reduction. Indicates rounding up;

[0059] Let T be the noise threshold of the edge map of the digital surface model. DSM The following formula is used to suppress noise in the edge map of the digital surface model:

[0060]

[0061] in, This is the edge map of the digital surface model.

[0062] Furthermore, the process of registering the digital surface model with the hyperspectral image using image registration methods includes the following steps:

[0063] The edge map of the denoised hyperspectral image in step 3. As a moving image, the edge map of the digital surface model As a fixed image, a feature extraction-based image registration method is used to register the two images, generating a two-dimensional affine transformation relationship T. affine ;

[0064] Let the hyperspectral images be [HSI1(x,y),HSI2(x,y),...,HSI S [x,y], image resolution is X×Y, according to Taffine Establish the pixel correspondence between the X×Y hyperspectral image and the digital surface model R×C image. Registration has been achieved.

[0065] Furthermore, after determining the grayscale of the raster in S2.2, median filtering needs to be performed on all rasters.

[0066] Beneficial effects:

[0067] (1) This invention proposes a point cloud position correction method based on ray tracing. First, by establishing a scanning ray model for pushbroom hyperspectral imaging, the coverage area of ​​a single hyperspectral scanning frame is determined, thereby obtaining the correspondence between the hyperspectral data frame and the lidar point cloud. Based on this principle, a fast point cloud layering algorithm is proposed to reduce its time complexity. Simultaneously, a point cloud position correction method based on ray tracing is proposed, generating an offset point cloud that can be used to generate a digital surface model containing stereoscopic information through orthographic projection, thus avoiding the occlusion detection process of curved surface projection.

[0068] (2) This invention also proposes a hyperspectral image edge feature map denoising algorithm. Based on the gray-scale distribution characteristics of feature points of hyperspectral single-band edge maps, a full-band cyclic denoising method is proposed to make up for the shortcomings of low spatial resolution of hyperspectral images by using the rich spectral dimensions of hyperspectral images.

[0069] Based on the features of this invention, it can effectively solve the problem of significant loss of stereoscopic dimensional information due to low-altitude viewing angle errors in stereo registration of hyperspectral images and lidar point clouds on low-altitude airborne platforms. Applying this invention to stereo registration of hyperspectral images and lidar point clouds on low-altitude airborne platforms can effectively improve the accuracy of registration and has a promising application prospect on low-altitude airborne platforms. Attached Figure Description

[0070] Figure 1 This is a flowchart illustrating the overall process of stereo registration of pushbroom hyperspectral image-LiDAR data based on ray tracing correction.

[0071] Figure 2 This is a flowchart of a fast algorithm for LiDAR point cloud correction based on ray projection.

[0072] Figure 3 This is a flowchart of a full-band cyclic noise reduction algorithm for hyperspectral image edge maps.

[0073] Figure 4 This is an example of attitude stability condition (a) for straight flight.

[0074] Figure 5 This is an example of attitude stability condition (b) for straight flight.

[0075] Figure 6 This is a layered diagram of the point cloud effect.

[0076] Figure 7 This is a diagram showing the effect of point cloud offset based on ray tracing.

[0077] Figure 8 A comparison of the generated digital surface model with the average hyperspectral image across the entire band.

[0078] Figure 9 A comparison of the edge map of the generated digital surface model and the edge map of the hyperspectral full-band model.

[0079] Figure 10 The image shows the result before registration.

[0080] Figure 11 This is a picture showing the result after registration.

[0081] Figure 12 The image shows the result after registration (the intersection-union area is marked).

[0082] Figure 13 This is a stereo registration effect diagram of the present invention.

[0083] Figure 14 This is a stereo registration effect diagram of the present invention.

[0084] Figure 15 This is a partial effect diagram of the stereo registration of the present invention.

[0085] Figure 16 This is a non-stereo registration effect diagram.

[0086] Figure 17 This is a partial effect image of non-stereo registration. Detailed Implementation

[0087] Specific implementation method one: Combining Figure 1 This implementation method is described below.

[0088] This embodiment is a pushbroom-based hyperspectral image-LiDAR data stereo registration method based on ray tracing correction, which includes the following process:

[0089] Step 1: Layer the lidar point cloud and establish the correspondence between lidar points and hyperspectral frames; for example... Figure 2 As shown, the process of layering lidar point clouds includes the following steps:

[0090] S1.1 To match the lidar point cloud with a single-frame hyperspectral image, it is necessary to calculate the layered planar information of the corresponding single frame using airborne positioning and attitude data. First, the spatial position and flight direction vector of the UAV during the scanning of the hyperspectral single frame need to be calculated:

[0091] Let the trajectory data be {P}m ,T m POS m}, where m = 1, 2, ..., M, M represents the number of hyperspectral image frames; P m =(E m N m H m () represents the coordinates of the E, N, H axes in the UTM coordinate system; T m For timestamps; ω is the direction angle of flight. m ,κ m , These are pitch angle, heading angle, and yaw angle, respectively.

[0092] Let d be the normalized vector of the flight direction at the scan time. m =A m ·v E +B m ·v N +C m ·v H , where v E ,v N ,v H Represents the unit direction vectors of the E, N, and H axes in UTM coordinates; A m B m C m Let be the normalization coefficients, and let the normalization coefficients satisfy .

[0093] Based on the flight heading and pitch angles at that moment, the direction vector at the frame scan moment is calculated using the following formula:

[0094]

[0095] S1.2, Spatial layer Z corresponding to a single hyperspectral scan frame m Represented as front and rear layered planes Π m Π m+1 Intermediate mezzanine space, single-layered plane Π i From the plane normal vector d i With a fixed point P on the plane i Determine, where i = 1, 2, ..., M, M+1.

[0096] Among them, planes Π2 to Π M The corresponding plane normal vector is the flight direction vector d at the previous and next time steps. m The average vector, with the fixed point being the flight coordinates P at the previous and next time moments. m The midpoint. The first and last planes Π1 and Π M+1 The normal vector is taken as the flight direction vectors d1 and d2 at the beginning and end of the flight time. MThe fixed points on the plane are Π2, Π3 and Π M-1 Π M The fixed point is extended in the opposite direction to find a point of equal length to determine the point.

[0097] Based on the determined plane normal vector d i Using the plane normal vector d i =A i ·v E +B i ·v N +C i ·v H Obtain coefficient A i B i C i Based on a fixed point P i =(E i N i H i ), (E i N i H i The coordinates of the E, N, and H axes in the UTM coordinate system, and the plane Π i as follows:

[0098] Π i :π i =A i ·e+B i ·n+C i ·h-(E i +N i +H i ) = 0

[0099] Where, π i =A i ·e+B i ·n+C i ·h-(E i +N i +H i ) represents the plane Π i The corresponding function, π i (p) represents substituting the coordinates of a point p = (E, N, H) in space into the function π. i The calculated function values; e, n, and h are the coordinates of the E, N, and H axes in the three UTM coordinate system, expressed as independent variables in the function form.

[0100] Based on M+1 boundary planes, the point cloud space is divided into M scanning layers Z corresponding to the hyperspectral image frames. m Each scan layer consists of a plane Π m With Π m+1 Determine any point p in the point cloud. n =(E n N nH n In the scanning layer Z m The necessary and sufficient condition for the inner can be expressed as:

[0101]

[0102] in, This indicates that the two are equivalent and are both necessary and sufficient conditions for each other.

[0103] Based on the order of frame scanning, the corresponding formulas for each point in the spatial layer are substituted for verification. For points p that satisfy the formulas... n The corresponding scanning layer number is recorded, but it is not used for verification of subsequent layers. Based on the layer plane information, the correspondence between each LiDAR point cloud and the spatial layer is calculated.

[0104] When the UAV's flight satisfies the stability conditions (a) and (b) for straight-line flight, a fast point cloud layering algorithm can be applied. That is, the processing in step 1.2 above can be replaced by a fast point cloud layering algorithm, thereby simplifying the process and reducing the computational load. The conditions for the UAV's flight to satisfy the stability conditions for straight-line flight and the fast point cloud layering algorithm are as follows:

[0105] (a) The flight trajectory can approximate a straight path, that is, it satisfies the instantaneous heading angle κ of all trajectory points in the path. m It tends to a constant K, with a deviation of no more than ±5°.

[0106] (b) During a single-segment straight flight, the altitude fluctuation does not exceed 5 meters, and the flight pitch angle ω is constant at all scanning times. m They approach 0 and none exceed 10°.

[0107] When the flight path is relatively stable and approximately straight, at this point, the layered planes are perpendicular to the ground surface, and all planes are parallel to each other, i.e., C i =0, A i B i Let A and B be constants, then plane Π i It can be obtained from a fixed point P i Confirmed, indicated as:

[0108] Π i :π i =A·e+B·nD i =0

[0109] The constant term D of the plane equation i =E i +N i +H i Let E, N, and H be the coordinates of fixed points, and let A = sinK and B = cosK.

[0110] During straight-line flight, the plane group constant term D iIt is a monotonic sequence. It passes through any point p in space. n And a plane parallel to the plane group can be represented as vtp n :A·e+B·nD n =0. Where the constant term D n Let E, N, and H be the coordinates of this point.

[0111] At this time p n In scan layer Z m The necessary and sufficient condition for the inner can be expressed as:

[0112]

[0113] Where D m D m+1 For scanning layer Z m The corresponding front and rear layer plane equations Π m Π m+1 The constant term.

[0114] When the hyperspectral imager flight data meets the straight-line attitude stability condition, the above method is used to first traverse the point cloud to obtain preliminary point cloud layering results p. n ∈Z m Then, verify the correctness of the stratification result using the non-approximate plane equation in equation (2). If the verification is incorrect, then p n Substitute the nearest neighbor layer Z m±t Verify according to formula (2) until the verification is correct.

[0115] Since this algorithm does not require traversing the point cloud in the calculation of each layer (in practical applications, the total number of points in the point cloud is much larger than the total number of layers), this method can reduce the time complexity of the algorithm from O(N^2) to O(N^2). 2 ) decreased to o(N).

[0116] Step 2: Based on the layered results, generate a digital surface model based on the offset projection of the ray-traced point cloud and establish its correspondence with the original point cloud. The process of generating the digital surface model based on the offset projection of the ray-traced point cloud includes the following steps:

[0117] S2.1 Point cloud offset based on ray tracing:

[0118] Based on the point cloud layering results from step 1, and based on any point p within the point cloud layer corresponding to a certain frame of the hyperspectral image... n =(E n N n H n ) Obtain the corresponding offset point Calculate using the following formula:

[0119]

[0120] Where κ' is the track angle, representing the angle relative to true north. N The acute angle formed; H pj P represents the height of the hyperspectral imaging plane. m =(E m N m H m () represents the spatial coordinates of the UAV at the moment of scanning this frame.

[0121] Iterate through all point clouds within the layer, perform offset transformations on these point clouds, and preserve their height H. n Without changing, generate an offset point cloud.

[0122] This step aims to reflect LiDAR points that match the stereo information of hyperspectral images in a digital surface model by offsetting the point cloud, that is, to offset and correct the point cloud by simulating the imaging angle of hyperspectral imaging.

[0123] S2.2 Establishing a digital surface model (LiDAR point cloud rasterization):

[0124] Create a raster with R×C rows and columns, where R×C should be greater than the resolution of the hyperspectral image. Based on the coordinate raster parameters, establish the correspondence between each point in the offset point cloud and the raster using orthographic projection.

[0125] in Indicates lidar points The coordinates of the corresponding grid in the digital surface model; (r,c) represents the horizontal and vertical coordinates of a certain grid in the lidar digital surface model.

[0126] The value of a single grid cell is the highest among all the height values ​​of all points within the cell. After traversing all grid cells, the grayscale (height value) is linearly mapped to the range [0,1], thus obtaining the digital surface model.

[0127] If the point cloud density is insufficient or the grid parameters are too large, resulting in some grids having no corresponding LiDAR points, the generated digital surface model may have salt-and-pepper noise. The preferred solution is to use median filtering for noise reduction. In fact, if the point cloud density is sufficient or all grids have corresponding LiDAR points, median filtering can also be used for noise reduction.

[0128] Step 3: Based on the principle of phase consistency, perform edge extraction on the digital surface model and the hyperspectral image, such as... Figure 3 As shown, it includes the following steps:

[0129] Single-band images in hyperspectral images With digital surface model I DSMAll images are grayscale, where S represents the total number of bands in the hyperspectral image. Two-dimensional Log-Gabor filters with 30° intervals are used to filter these grayscale images. Edge maps M are generated using a phase consistency algorithm. These edge maps are categorized into groups, denoted as: The edge graph of the digital surface model is M DSM .

[0130] Let the noise threshold group of the highlight image edge group be . threshold Determined by the grayscale distribution of each edge image. Noise suppression is performed by traversing the entire band edge image of the specular image using the following formula:

[0131]

[0132] in, This is the edge map of the hyperspectral image after full-band noise reduction. This indicates rounding up to the nearest integer.

[0133] Let T be the noise threshold of the edge map of the digital surface model. DSM The following formula is used to suppress noise in the edge map of the digital surface model:

[0134]

[0135] This is the edge map of the digital surface model.

[0136] This step aims to improve the consistency between hyperspectral images and digital surface models through a unified edge extraction algorithm, thereby facilitating the application of image registration methods in step 4.

[0137] Step 4: Register the digital surface model with the hyperspectral image using image registration methods to generate a hyperspectral point cloud, including the following steps:

[0138] The edge map of the denoised hyperspectral image in step 3. As a moving image, the edge map of the digital surface model As a fixed image, a feature extraction-based image registration method is used to register the two images, generating a two-dimensional affine transformation relationship T. affine .

[0139] Let the hyperspectral images be [HSI1(x,y),HSI2(x,y),...,HSI S [x,y], image resolution is X×Y, according to T affine Establish the pixel correspondence between the X×Y hyperspectral image and the digital surface model R×C image. Registration has been achieved.

[0140] The process of this invention can be simply represented as the correspondence between the offset point cloud and the digital surface model in step 2. The hyperspectral image pixels (x, y) are compared with the original point cloud p. n Correspondence. The process of establishing the correspondence between lidar point clouds and hyperspectral image pixels is shown in the following formula:

[0141]

[0142] Example

[0143] To verify the effectiveness of the present invention, hyperspectral image-point cloud stereo registration based on ray tracing correction was performed using the process of Specific Embodiment 1.

[0144] Examples of the "straight flight attitude stability conditions" proposed in step 1 are as follows: Figures 4-5 As shown.

[0145] In this embodiment, the flight altitude fluctuates by about 1 meter, the flight path is approximately a straight line, the instantaneous heading angle deviation does not exceed ±5°, and the flight pitch angle approaches 0.

[0146] The point cloud layering results obtained in step 1 are shown in Figure 6.

[0147] The point cloud offset results based on ray tracing in step 2 are as follows: Figure 7 As shown, gray (actually red in the color image) represents the original point cloud, and black (actually blue in the color image) represents the corrected point cloud. It can be observed that the corrected point cloud shows a shift away from both sides of the flight path, with the shift at higher altitudes being greater than that at lower altitudes.

[0148] A comparison of the digital surface model generated in step 2 (left) and the hyperspectral full-band mean image (right). Figure 8 As shown, the digital surface model contains some three-dimensional planar information, similar to the hyperspectral image, but the two image features differ significantly.

[0149] The edge map (left) and hyperspectral full-band edge map generated in step 3 are shown below. Figure 9 As shown, the consistency of features between the two can be observed to be significantly improved.

[0150] The before-and-after registration results between the hyperspectral image and the digital surface model in step 4 are shown below. Figures 10-12 As shown, Figure 10 To show the effect before registration, Figure 11 This is a picture showing the result after registration. Figure 12 This is the image after registration (with the intersection-over-union (IoU) region marked). It can be observed that the line feature matching degree and IoU of the key regions in the registered image are high.

[0151] The effect after stereo registration in step 4 is as follows: Figures 13-15 As shown, where Figure 13 This is a stereo registration effect diagram of the present invention. Figure 14 This is a stereo registration effect diagram of the present invention. Figure 15 This is a partial effect diagram of the stereo registration of the present invention.

[0152] The "stereo" registration of this invention is reflected in the point cloud correction based on ray tracing in step 2. To illustrate the "stereo" effect of this invention, a non-stereo registration is also performed on the scheme that removes step 2. The non-stereo registration effect is as follows: Figures 16-17 As shown, where, Figure 16 This is a non-stereo registration effect image. Figure 17 This is a partial effect image of non-stereo registration.

[0153] In non-stereo registration, it is impossible to register two types of data in a stereo dimension (such as building walls); however, the stereo registration method of this invention can correctly achieve the registration of stereo dimension information.

[0154] The above examples of the present invention are merely illustrative of the computational model and process of the present invention, and are not intended to limit the implementation of the present invention. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is impossible to exhaustively list all possible implementations here. Any obvious variations or modifications derived from the technical solutions of the present invention are still within the scope of protection of the present invention.

Claims

1. A hyperspectral image-point cloud stereo registration method based on ray tracing correction, characterized in that, Includes the following steps: S1. Layer the lidar point cloud, including the following steps: S1.1 Calculate the spatial position and flight direction vector of the UAV when scanning a single hyperspectral frame: Set trajectory data ,in M represents the number of hyperspectral image frames; Represents the coordinates of the E, N, and H axes in the UTM coordinate system; For timestamps; The direction angle of flight, These are pitch angle, heading angle, and yaw angle, respectively. The direction vector at each frame scan time is calculated using the following formula: （1） in, Represents the unit direction vectors of the E, N, and H axes in UTM coordinates; S1.2, Spatial layer of lidar point cloud corresponding to a single hyperspectral scanning frame Represented as front and rear layered planes , Intermediate mezzanine space, based on layered plane Determine the spatial layer of the lidar point cloud corresponding to a single hyperspectral scan frame. ; Layered plane From the plane normal vector With a fixed point on the plane Confirmed, among which ; S2. Based on the layered results, generate a digital surface model based on the offset projection of ray-traced point clouds and establish a correspondence with the original point cloud. The process of generating a digital surface model based on the offset projection of ray-traced point clouds includes the following steps: S2.1 Based on the point cloud layering results from step 1, and based on any point within the point cloud layer corresponding to a certain frame of the hyperspectral image... Obtain the corresponding offset point Calculated according to the following formula: in, The track angle indicates the direction relative to true north. The acute angle formed; The height of the hyperspectral imaging plane; To scan the spatial coordinates of the UAV at that moment in the frame; Iterate through all point clouds within the layer, perform offset transformations on these point clouds, and preserve their height. Without changing anything, generate an offset point cloud; S2.2, Establish a system with the following number of rows and columns. The grid, The resolution should be greater than that of the hyperspectral image; the offset point cloud should be orthographically projected to establish the correspondence between each point and the grid based on the coordinate grid parameters. ;in, Indicates lidar points The coordinates of the corresponding grid in the digital surface model; Represents the horizontal and vertical coordinates of a grid cell in the lidar digital surface model; The value of a single raster cell is the highest height value among all points within that cell, denoted as the raster's grayscale. After traversing all raster cells, the grayscale is linearly mapped to... Inside, a digital surface model is obtained; S3. Based on the principle of phase consistency, edge extraction is performed on the digital surface model and the hyperspectral image; S4. Using image registration methods, the digital surface model is registered with the hyperspectral image to generate a hyperspectral point cloud.

2. The hyperspectral image-point cloud stereo registration method based on ray tracing correction according to claim 1, characterized in that, Layered plane in S1.2 normal vector on the plane Determined in the following manner: Layered plane to The corresponding plane normal vector is the flight direction vector at the previous and next time moments. The average vector; First and last plane and The normal vector is taken as the flight direction vector at the beginning and end of the time. and .

3. The hyperspectral image-point cloud stereo registration method based on ray tracing correction according to claim 2, characterized in that, Layered plane in S1.2 Fixed points on Determined in the following manner: Layered plane to The fixed point on the plane is the flight coordinate of the previous and next time moments. The midpoint; First plane The plane fixed point is formed by The fixed point is extended in the opposite direction to find a point of equal length to determine the position; end plane The fixed point of the surface is determined by The fixed point is extended in the opposite direction to find a point of equal length to determine the point.

4. The hyperspectral image-point cloud stereo registration method based on ray tracing correction according to claim 3, characterized in that, The layered plane described in S1.2 From the plane normal vector With a fixed point on the plane The determination process includes the following steps: Based on the determined plane normal vector Using plane normal vectors Obtain the coefficient Based on fixed points , The coordinates of the E, N, and H axes in the UTM coordinate system, and the plane as follows: in, Representing a plane The corresponding function, This indicates that a point in space Substitute the coordinates into the function The calculated function value; The coordinates of the E, N, and H axes in the three UTM coordinate system are expressed as independent variables in the form of a function.

5. The hyperspectral image-point cloud stereo registration method based on ray tracing correction according to claim 4, characterized in that, Based on layered plane Determine the spatial layer corresponding to a single hyperspectral scan frame The process includes the following steps: according to A dividing plane divides the point cloud space into M scanning layers corresponding to the hyperspectral image frames. Each scan layer consists of a plane and Determine any point in the point cloud. In the scanning layer The necessary and sufficient condition within is expressed as: （2） in, This indicates that the two are equivalent; Based on the order of frame scanning, each point is substituted into the corresponding formula of the spatial layer for verification. For points that satisfy the formula... Record the corresponding scanning layer number, which is not used for subsequent layer verification; calculate the correspondence between each lidar point cloud and the spatial layer based on the layer plane information.

6. The hyperspectral image-point cloud stereo registration method based on ray tracing correction according to claim 3, characterized in that, In step 1.2, if the UAV flight satisfies the straight-line flight attitude stability conditions (a) and (b), the fast point cloud layering algorithm is used to determine the layering plane. and based on the hierarchical plane Determine the spatial layer corresponding to a single hyperspectral scan frame ; The drone's flight meets the following attitude stability conditions for straight-line flight: (a) The flight trajectory can approximate a straight path, that is, it satisfies the instantaneous heading angle of all trajectory points in the path. tending to a constant And the deviation does not exceed ±5°; (b) During a single straight flight segment, the altitude fluctuation does not exceed 5 meters, and the flight pitch angle at all scanning moments is... Approaching 0, and all not exceeding 10°; Determine the layering plane using a fast point cloud layering algorithm The process includes the following steps: When the flight path is relatively stable and approximately straight, at this point, the layered planes are perpendicular to the ground surface, and all planes are parallel to each other. , Let A and B be constants, then the plane at this time From a fixed point Confirmed, indicated as: Among them, the constant term of the plane equation Let E be the coordinates of the fixed point E, N, and H axes. , .

7. The hyperspectral image-point cloud stereo registration method based on ray tracing correction according to claim 6, characterized in that, Based on layered plane Determine the spatial layer corresponding to a single hyperspectral scan frame Includes the following steps: In straight-line flight, the constant term of the plane group It is a monotonic sequence; it passes through any point in space. And the plane parallel to the plane group is represented as ; where the constant term Let E, N, and H be the coordinates of this point. at this time In the scanning layer The necessary and sufficient condition within is expressed as: （3） in, For scanning layer Corresponding front and back layer plane equations , constant term; When the hyperspectral imager flight data meets the straight-line attitude stability condition, the point cloud is first traversed to obtain preliminary point cloud layering results. And then according to The non-approximate plane equations in the model are used to verify the correctness of the stratification results; if the verification is incorrect, then... Bring in the nearest neighbor layer Continue until it is verified to be correct.

8. The hyperspectral image-point cloud stereo registration method based on ray tracing correction according to any one of claims 1 to 7, characterized in that, The process of edge extraction from digital surface models and hyperspectral images includes the following steps: Single-band images in hyperspectral images With digital surface model All images are grayscale images, where S represents the total number of bands in the hyperspectral image. A two-dimensional Log-Gabor filter is used to filter these grayscale images, and an edge map M is generated using a phase consistency algorithm. The edge map group of the hyperspectral image is denoted as... The edge map of the digital surface model is ; Let the noise threshold group of the highlight image edge group be . Noise suppression is achieved by iterating through the entire band edge map of the specular image using the following formula: （5） in, This is the edge map of the hyperspectral image after full-band noise reduction. Indicates rounding up; Let the noise threshold of the edge map of the digital surface model be . The following formula is used to suppress noise in the edge map of the digital surface model: （6） in, This is the edge map of the digital surface model.

9. The hyperspectral image-point cloud stereo registration method based on ray tracing correction according to claim 8, characterized in that, The process of registering a digital surface model with a hyperspectral image using image registration methods includes the following steps: The edge map of the denoised hyperspectral image in step 3. As a moving image, the edge map of the digital surface model As a fixed image, a feature extraction-based image registration method is used to register the two images, generating a two-dimensional affine transformation relationship. ; The hyperspectral image is denoted as The image resolution is ,according to Establish Hyperspectral image to digital surface model Pixel correspondence of the image This means that registration has been achieved.

10. The hyperspectral image-point cloud stereo registration method based on ray tracing correction according to claim 9, characterized in that, After determining the grayscale of the raster in S2.2, median filtering needs to be performed on all rasters.

Images