Three-dimensional multispectral point cloud reconstruction method, system and device based on dual-angle multispectral image

Abstract

The application discloses a three-dimensional multi-spectrum point cloud reconstruction method, system and device based on a double-angle multi-spectrum image, and belongs to the technical field of multi-spectrum image three-dimensional reconstruction. In order to solve the problem that the existing three-dimensional reconstruction method cannot directly reconstruct complete waveband point clouds from non-registered cameras, the application first performs texture enhancement on images collected by cameras with different angles; the cameras are multi-spectrum cameras; tilt image projection correction and feature point extraction are performed on the tilt images, and then the feature point coordinates are inversely transformed back to the tilt images, so that the feature points in the tilt images are determined, and the original images are matched and searched for image pairs and image matching based on the feature points; then, multi-view dense matching is performed on the images to obtain a depth map point cloud and a multi-spectrum depth map point cloud, the multi-spectrum depth map point clouds collected by the same camera are spliced, the point clouds corresponding to different multi-spectrum cameras and containing different waveband spectrums are spliced to obtain full-waveband multi-spectrum point clouds, and finally, abnormal points are filtered out.

Description

Technical Field

[0001] This invention belongs to the field of multispectral image three-dimensional reconstruction technology, specifically relating to a three-dimensional multispectral point cloud reconstruction method, system, storage medium and device based on dual-angle multispectral images. Background Technology

[0002] With the development of sensor and drone technologies, drone remote sensing has gained widespread application in many fields, including military defense, urban management, land and environmental monitoring, and agricultural informatization, due to its advantages such as high mobility, high real-time performance, and avoidance of cloud cover. Drone platforms can be equipped with different sensors according to specific mission requirements. Among them, multispectral cameras can acquire multi-band spectral images with high spatial resolution, while also offering advantages such as low cost and high efficiency, playing a vital role in applications such as agriculture, ecology, and urban monitoring.

[0003] Limited by the two-dimensional spatial dimension of images, remote sensing observation methods based on multispectral images cannot accurately and completely describe the spatial geometric features of the observed scene. Specifically, they cannot perceive the true size and height of ground objects, leading to a lack of accuracy in many quantitative remote sensing analysis applications. 3D reconstruction is a technique that uses multiple overlapping images to recover the three-dimensional spatial information of an observed scene. It can eliminate geometric distortions and ground object occlusion in images, allowing for analysis of the observed scene from a three-dimensional spatial perspective. Point clouds are the primary target data for 3D reconstruction algorithms. They are collections of a large number of data points, each of which possesses three-dimensional spatial information and can incorporate multidimensional spectral information, effectively compensating for the deficiencies of image data.

[0004] Vertical ground observation is currently the primary method for acquiring multispectral images from UAVs, but it lacks the ability to obtain side-view images of ground features in the observed scene, resulting in low integrity of 3D reconstruction results. With the diversification of sensor technology, equipping multiple tilting cameras to capture images of the ground from different angles can alleviate the problem of insufficient acquisition of side-view images of ground features and play a significant role in 3D reconstruction tasks.

[0005] Based on the above problems and research background, this invention combines the advantages of multispectral imaging technology, three-dimensional reconstruction technology and oblique photography technology. It uses a dual-angle multispectral camera to acquire two sets of five-band multispectral images and reconstructs three-dimensional multispectral point clouds from the dual-angle multispectral images. This has important research significance and application value for qualitative and quantitative remote sensing such as spectral index analysis, three-dimensional scene modeling, and multi-dimensional stereo detection. Summary of the Invention

[0006] This invention aims to address the problem that existing 3D reconstruction methods cannot directly reconstruct complete band point clouds from unregistered cameras.

[0007] A method for reconstructing 3D multispectral point clouds based on dual-angle multispectral images includes the following steps:

[0008] S1. Perform texture enhancement on images acquired by cameras with different set angles; the camera is a multispectral camera, and the spectral bands corresponding to cameras with different set angles are different, and the images acquired by the multispectral camera are multispectral reflectance maps;

[0009] S2. Based on the texture-enhanced image of each multispectral image, for the tilted image, tilted image projection correction is performed. The scale-invariant feature transform matching algorithm and SIFT algorithm are used to extract feature points from the projection-corrected image. Then, the coordinates of the feature points are inversely transformed back to the tilted image to determine the feature points in the tilted image, which are also the feature points in the original image.

[0010] S3. Based on the feature points in the original image, perform a matching image pair search on the original image, and search for images whose overlap rate with the reference image exceeds the overlap rate threshold for subsequent matching.

[0011] S4. Perform image matching on the images found in S3;

[0012] S5. Based on the matching point pairs determined in S4, estimate the extrinsic parameters of each image and the intrinsic parameters of the camera, and perform multi-view dense matching on the images to obtain the depth map point cloud.

[0013] After obtaining the depth map point cloud for each image, the multispectral reflectance map is distorted by combining the camera intrinsic parameters and image extrinsic parameters, and the multiband spectral values ​​of each point in the depth map point cloud are obtained from it, thus obtaining the multispectral depth map point cloud.

[0014] Then, the multispectral depth point clouds collected by the same camera are stitched together to obtain multispectral depth point clouds of the same band. Then, based on the consistency of spatial position, the point clouds of different band spectra corresponding to cameras with different settings angles are stitched together to generate a complete full-band multispectral point cloud. Finally, outliers are filtered out.

[0015] The aforementioned full-band multispectral point cloud is a multispectral point cloud covering all bands of different cameras.

[0016] Furthermore, texture enhancement for images acquired by different multispectral cameras includes the following steps:

[0017] First, images are acquired using different multispectral cameras, and then the Normalized Differential Vegetation Index (NDVI) is calculated using each camera.

[0018] Define the operation rule A = [B ⊙ λ b ] indicates that for each element in matrix B, if its value ⊙λ bIf λ = 1, then the corresponding element in matrix A with the same dimension as B is 1; otherwise, it is 0. b ⊙ represents a numerical value, and ⊙ represents one of the symbols >, ≥, <, or ≤.

[0019] The following formula can be used to remove the shaded areas from the vegetation index map.

[0020]

[0021] I NDVI NDVI, a vegetation index map used to represent different multispectral cameras. j A single-band reflectance map used to represent the center wavelength j of different multispectral cameras;

[0022] Then, the vegetation pixels in the vegetation index map are combined with the non-vegetation pixels in the blue, green, and red band reflectance maps. Pixels with a vegetation index value less than 0.3 are considered non-vegetation pixels. The calculation formula is as follows:

[0023]

[0024] In the formula, Image I represents the reflectance images of the corresponding red, green, and blue bands in the multispectral image, respectively; Image I is the final NDVI-RGB composite image with strong geometric texture.

[0025] Furthermore, the process of performing tilted image projection correction includes the following steps:

[0026] First, determine the perspective projection transformation relationship corresponding to the downward-viewing image after correcting the tilted image. Let p(x,y,1) and P′(X,Y,Z) represent the image coordinates before perspective projection and the spatial coordinates after projection, respectively.

[0027]

[0028] Where M is the image perspective projection matrix, m 11 to m 33 For elements in the image perspective projection matrix;

[0029] Let p′(x′,y′) be the coordinates of the projected image. From equation (3), we can establish the equation shown in equation (4).

[0030]

[0031] Let m 33 =1, then the number of unknown variables in the equation becomes 8, and matrix M can be solved by using 4 sets of image coordinates before and after projection; using the coordinates of the 4 vertices of the tilted image as p, the corresponding 4 projection points p′ are solved:

[0032] First, project the homogeneous pixel coordinates o(u0,v0,1) of the center point of the tilted image backward onto the horizontal ground. Let the world coordinates of its projection point be... The rotation matrix R1 of the tilted image in the world coordinate system can be obtained from the pose data.

[0033] Then, the latitude, longitude, and altitude coordinates of each camera exposure position are converted into geocentric fixed rectangular coordinates to obtain the translation vector t1 of the image in the world coordinate system;

[0034] Based on the rotation matrix R1 of the image in the world coordinate system, the ground projection point of the image center point is obtained according to equation (7):

[0035]

[0036] In the formula, K1 represents the camera intrinsic parameters of the tilted image; D represents the height from the camera optical center to the ground, i.e., the depth of the projection point; let the initial depth D be the flight altitude of the UAV, when... Time-updated depth And recalculate using equation (7), where α1 is the world coordinate Z-axis reference threshold, until... The coordinates of the ground projection point can be obtained at that time. Simultaneously obtain its depth D;

[0037] The ground point coordinates are obtained by projecting the vertex coordinates of the tilted image onto the ground.

[0038] Then, set the roll and pitch angles of the corrected image to 0, and keep the yaw angle the same as that of the tilt image, i.e., (α,β,γ)=(0,0,γ). corr ), and obtain the rotation matrix R2 of the corrected image in the world coordinate system;

[0039] Let the corrected image be at world coordinates o. C (X W ,Y W Let h = D, then we get the translation vector t2 of the corrected image;

[0040]

[0041] Finally, use equation (9) to combine the four ground points. By projecting the image forward onto the corrected image, we can obtain the corresponding four image projection points.

[0042] p′=K1[R2p g +t2] (9)

[0043] Substituting the four sets of image coordinates before and after projection into equation (4) yields eight equations. Solving these equations provides the image projection transformation matrix M, which is then applied to equation (3) to transform the tilted image. The transformed spatial coordinates P′(X,Y,Z) are then normalized to image coordinates, i.e. This completes the projection correction of the tilted image.

[0044] Furthermore, the process of performing image pair matching search on the original image based on feature points in the original image includes the following steps:

[0045] First, the geocentric fixed rectangular coordinates of each original image are obtained based on its pose data. Simultaneously determine the minimum longitude (lon) in the original image within this dataset. ref latitude ref Let the reference height be alt. ref Let the point in space where 0 = 0 be the reference point, and find the geocentric fixed rectangular coordinates of the reference point. The coordinates of each original image in the world coordinate system are obtained from equation (10), and the rotation matrix R and translation vector t in the world coordinate system are solved respectively.

[0046]

[0047] Then, the center point of each image is projected backwards onto the ground to obtain the projection point O. G (X G ,Y G ,0);

[0048] Suppose a certain camera is at O C The ground projection coordinates of the center point of a rear-view oblique image taken at the location are O. G Then the neighbor matching image projection point O′ of the reference image G Satisfying equation (11):

[0049] ||O G -O′ G ||≤r s (11)

[0050] Where, r s The search radius is estimated based on the image overlap rate.

[0051] Images whose overlap with the reference image exceeds an overlap threshold are searched for for subsequent matching.

[0052] Furthermore, the process of image matching for the images searched by S3 includes the following steps:

[0053] First, the KD-Tree algorithm is used to search for the nearest and second nearest neighbors of the feature points in the reference image on the matching image. Let the Euclidean distances between the nearest and second nearest neighbors and the feature points in the reference image be d1 and d2, respectively. When the two satisfy equation (13), the feature points in the reference image match the nearest neighbors, that is, r is used as the confidence of the matching point pair.

[0054]

[0055] Where η is the matching point filtering threshold;

[0056] For the matching results, the MAGSAC++ algorithm is used to remove mismatched point pairs. In the process of removing mismatched point pairs, the homography transformation matrix is ​​estimated by maximizing the consistency principle.

[0057] The homography transformation matrix M between the two images to be matched is estimated. d Next, the coordinates of the feature points in the matching image are projected onto the coordinate system of the reference image, and the reference image is divided into blocks, each block being a sub-image. Then, each sub-image of the reference image is traversed, and local re-matching is performed, that is, KD-Tree nearest neighbor search matching is performed on the reference feature points and matching feature points within the sub-image region, and then the matching results of all sub-images are merged.

[0058] Swap the reference image and the matching image, while utilizing M d -1 Project the feature points of the new matching image, divide the new reference image into blocks, and perform local rematching.

[0059] The intersection of the bidirectional matching results is taken as the matching result of the two images;

[0060] Finally, the MAGSAC++ algorithm is used to remove mismatched point pairs from the matching results to obtain the final matching point pairs.

[0061] Furthermore, in S5, the process of estimating the extrinsic parameters of each image and the intrinsic parameters of the camera based on the matching point pairs determined in S4 is achieved through the principle of motion reconstruction structure.

[0062] Furthermore, in the final outlier filtering process described in S5, outlier identification is performed using a statistical method from the open-source point cloud data processing library PDAL:

[0063] For each point in the multispectral point cloud, search for k nearest neighbors and calculate the average distance. The outlier detection threshold is obtained by summing the standard deviation σ.

[0064]

[0065] Where m is the multiplication factor;

[0066] When the distance u of the nearest neighbor of i i Points ≥ t will be marked as outliers and filtered out.

[0067] A 3D multispectral point cloud reconstruction system based on dual-angle multispectral images includes an image texture enhancement unit, a projection correction and feature point extraction unit, a matching image pair search unit, an image matching unit, and a full-band multispectral point cloud generation unit; among which...

[0068] Image texture enhancement unit: performs texture enhancement on images acquired by different multispectral cameras; the different multispectral cameras are cameras with different spectral bands, and the images acquired by the multispectral cameras are multispectral reflectance maps;

[0069] Projection correction and feature point extraction unit: Based on the texture-enhanced image of each multispectral image, for the tilted image, tilted image projection correction is performed, the SIFT algorithm is used to extract feature points from the projection-corrected image, and then the coordinates of the feature points are inversely transformed back to the tilted image, thereby determining the feature points in the tilted image, which are also the feature points in the original image;

[0070] Matching Image Pair Search Unit: Based on feature points in the original image, the original image is searched for matching image pairs. Neighboring images with an overlap rate exceeding the overlap rate threshold with the reference image are searched for and used for subsequent matching.

[0071] Image matching unit: performs image matching on the searched images;

[0072] Full-band multispectral point cloud generation unit: First, based on the determined matching point pairs, estimate the extrinsic parameters of each image and the intrinsic parameters of the camera, and perform multi-view dense matching on the images to obtain depth map point clouds;

[0073] After obtaining the depth map point cloud for each image, the multispectral reflectance map is distorted by combining the camera intrinsic parameters and image extrinsic parameters, and the multiband spectral values ​​of each point in the depth map point cloud are obtained from it, thus obtaining the multispectral depth map point cloud.

[0074] Then, the multispectral depth map point clouds collected by the same camera are stitched together to obtain multispectral depth map point clouds of the same band. Then, the point clouds containing different bands of spectrum from different multispectral cameras are stitched together to generate a complete full-band multispectral point cloud based on spatial location consistency. Finally, outliers are filtered out.

[0075] The aforementioned full-band multispectral point cloud is a multispectral point cloud covering all bands of different multispectral cameras.

[0076] A computer storage medium storing at least one instruction, which is loaded and executed by a processor to implement the aforementioned method for reconstructing a three-dimensional multispectral point cloud based on a dual-angle multispectral image.

[0077] A three-dimensional multispectral point cloud reconstruction device based on dual-angle multispectral images is disclosed. The device includes a processor and a memory. The memory stores at least one instruction, which is loaded and executed by the processor to implement the three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images.

[0078] Beneficial effects:

[0079] 1. The multispectral image texture enhancement method in this invention can improve the contrast between vegetation and non-vegetation pixels in an image, enhance image grayscale texture, and reduce multi-band grayscale nonlinear variation. The tilted image projection correction method can correct large-angle images to a downward viewing angle, significantly reducing the degree of perspective geometric distortion. The matching image pair search and local image re-matching method can significantly reduce mismatches and increase the number of correct matches. The above processing algorithms greatly improve the feature point extraction and matching logarithm of tilted multispectral images, ensuring the accuracy and integrity of subsequent 3D reconstruction algorithms.

[0080] 2. This invention uses a tilted dual-angle viewpoint to acquire multispectral images. Compared with a single-angle downward viewpoint, it can acquire more object side facade images. Combined with projection correction and local rematching methods, it can obtain more object side facade reconstruction point clouds, and the spatial structure of the reconstructed scene has high integrity.

[0081] 3. Since images captured by cameras from different perspectives cannot be directly registered, existing 3D reconstruction methods cannot directly reconstruct complete spectral point clouds from unregistered cameras. The depth map-based multispectral point cloud generation method designed in this invention first estimates the depth map point cloud of each image and merges the point clouds of the same camera. Then, it extracts the nearest neighbor points from all camera point clouds to obtain the complete spectral band reflectance value of each spatial point, thereby realizing the generation of complete spectral multispectral point clouds. Attached Figure Description

[0082] Figure 1 This is a basic flowchart of a three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images.

[0083] Figure 2 This is a schematic diagram showing that some wavebands of different ground features may have the same reflectivity.

[0084] Figure 3 This is a schematic diagram of the main processes of projection correction and feature extraction.

[0085] Figure 4 This is a schematic diagram illustrating the basic principle of projection correction.

[0086] Figure 5 A schematic diagram illustrating the basic principle of image pair search.

[0087] Figure 6(a) is a schematic diagram of mismatch in feature point matching; Figure 6(b) is a schematic diagram of missed match in feature point matching.

[0088] Figure 7 The flowchart shows the feature point matching process.

[0089] Figure 8 This is a schematic diagram of local feature point rematching.

[0090] Figure 9 This is a schematic diagram of the motion recovery structure algorithm.

[0091] Figure 10 This is a schematic diagram of the multispectral point cloud generation process.

[0092] Figure 11 To generate 10-band point clouds based on spatial location consistency.

[0093] Figure 12 This is a partial image capture.

[0094] Figure 13 This is a schematic diagram showing the texture enhancement results and a magnified comparison of a local area.

[0095] Figure 14 This is a comparison image before and after image projection correction.

[0096] Figure 15 A visualization of the search results for matching image pairs.

[0097] Figure 16 A visualization of a set of matching results.

[0098] Figure 17 This provides visualization of multispectral dense point clouds and results of vegetation index calculations.

[0099] Figure 18 This is the result of point cloud extraction for a single ground feature. Detailed Implementation

[0100] This invention is based on pose information-assisted dual-angle multispectral image feature extraction and matching, and then multispectral point cloud reconstruction based on depth map. It can reconstruct a three-dimensional point cloud model with complete spectral bands and fine spatial structure for dual-angle multispectral images of large scenes and complex land cover types.

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

[0102] This embodiment is a three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images. First, the experimental data used in this invention will be described. The main information is shown in Table 1. The pose sensor data is stored in the metadata area of ​​each image. The flight altitude is 80 meters, and the heading and flight strip overlap rate are both 80%.

[0103] Table 1 Main Information of Experimental Data

[0104]

[0105] The three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images described in this embodiment includes the following steps:

[0106] S1, Multispectral Image Texture Enhancement:

[0107] While multispectral images can reflect the spectral characteristics of ground features, the varying reflectance of ground features across different spectral bands leads to non-linear variations in pixel brightness for the same feature across different bands. Furthermore, single-band reflectance maps often exhibit weak texture, causing a sharp decline in the effectiveness of gradient-based feature extraction methods. For scenes containing various ground features such as buildings, vehicles, water bodies, trees, and bare land, the single-band reflectance map processing method is highly sensitive to band selection and changes in ground feature type. For example, two different ground feature points may have the same reflectance value in a certain band, making them indistinguishable during feature extraction and consequently affecting matching and reconstruction accuracy. Figure 2 As shown.

[0108] To address the aforementioned problems, this invention enhances image texture based on the spectral characteristics of ground features. First, images are acquired using two cameras, designated Camera 1 and Camera 2. In this embodiment, Camera 1 and Camera 2 are two tilted cameras (the cameras are tilted for the measured area). Then, the Normalized Differential Vegetation Index (NDVI) is calculated using the near-infrared or red-edge band and red band images from each camera. The calculation formula is as follows:

[0109]

[0110] In the formula, b i This is a single-band reflectance diagram with center wavelength i. These are vegetation index maps calculated from images of camera 1 and camera 2, respectively. Since camera 2 used in this invention does not have a near-infrared band, an image of the red-edge 740 band with a relatively close center wavelength is selected instead.

[0111] Define the operation rule A = [B > b] to mean that for each element in matrix B, if its value is greater than b, then the element at the corresponding position in matrix A of the same dimension as B has a value of 1, otherwise it has a value of 0; for other comparison relations similar to B > b, such as greater than or equal to, less than, less than or equal to, the same method is used.

[0112] The following formula is used to remove the shaded areas from the vegetation index map. Depending on the camera being processed, I NDVI b j They are respectively b 842 or b 740 .

[0113] I corr =I NDVI ×M shd

[0114] M shd =[b j ≥0.2]

[0115] Because vegetation in images typically has a lot of repetitive texture and its position is prone to change at different times, it is not conducive to image matching. Artificial features, on the other hand, have significant point features and fixed positions, making them the primary basis for image matching. However, in the vegetation index map, the index values ​​of non-vegetation pixels are small, and the image texture is not obvious. To further enhance the geometric texture and preserve the original texture features of artificial features, this invention synthesizes vegetation pixels in the vegetation index map with non-vegetation pixels in the blue, green, and red three-band reflectance maps. Pixels with a vegetation index value less than 0.3 are considered non-vegetation pixels. The calculation formula is as follows:

[0116]

[0117] In the formula Image I represents the reflectance images of the corresponding red, green, and blue bands in the multispectral image, respectively. Image I is the final NDVI-RGB composite image with strong geometric texture, and subsequent algorithms will extract and match feature points based on this composite image.

[0118] S2. Based on the texture-enhanced image of each multispectral image, perform tilt image projection correction for tilted images:

[0119] Although Camera 1 and Camera 2 are two tilted multispectral cameras, the data acquisition process is easily affected by the flight attitude of the UAV, which may result in the acquired images being flat. Even if Camera 1 or Camera 2 is not a tilted camera, the UAV's flight attitude may still result in the acquisition of tilted images. The tilted image projection correction process is for tilted images. If the image is not tilted, projection correction is optional.

[0120] like Figure 3 As shown, the process of correcting the tilted image to a downward-viewing image by performing projection correction and feature extraction is a perspective projection transformation, as shown in Equation (3), where p(x,y,1) and P′(X,Y,Z) represent the image coordinates before perspective projection and the spatial coordinates after projection, respectively, and M is the image perspective projection matrix.

[0121]

[0122] Where, m 11 to m 33 These are elements in the image perspective projection matrix.

[0123] Let p′(x′,y′) be the coordinates of the projected image (P′ is the coordinates of a point in three-dimensional space, and p′ is the coordinates of the corresponding two-dimensional image point). From equation (3), we can establish the equation shown in equation (4), and let m 33 If = 1, then the number of unknown variables in the equation becomes 8. Therefore, only 4 sets of image coordinates before and after projection are needed to solve matrix M.

[0124]

[0125] Therefore, using the coordinates of the four vertices of the tilted image as p, we solve for the corresponding four projection points p′:

[0126] First, based on the camera pose data and sensor parameters, the homogeneous pixel coordinates o(u0,v0,1) of the center point of the tilted image are projected backward onto the horizontal ground. Let the world coordinates of the projected point be... like Figure 4 As shown.

[0127] The rotation matrix R1 of the tilted image in the scene world coordinate system can be obtained using equation (5) based on the pose data:

[0128]

[0129] Where α, β, and γ are the roll angle, pitch angle, and yaw angle, respectively, and R Z (γ), R Y (β), R X (α) are two-dimensional rotation matrices about the Z, Y, and X axes, respectively;

[0130] Using equation (6), the latitude, longitude, and altitude coordinates (lon, lat, alt) of each camera exposure position are converted into geocentric fixed rectangular coordinates to obtain the translation vector t1(X, Y, Z) of the image in the world coordinate system of the scene:

[0131]

[0132] In the formula, N is the radius of curvature of the Earth, a = 6378137.0 is the major radius of the Earth's reference ellipsoid, and f = 1 / 298.257223565 is the extreme flattening of the Earth's reference ellipsoid.

[0133] Based on the rotation matrix R1 of the image in the scene world coordinate system, the ground projection point of the image center point is obtained according to equation (7):

[0134]

[0135] In the formula, D represents the height from the camera's optical center to the ground, i.e., the depth of the projection point. Let the initial depth D be the UAV's flight altitude, when... Time-updated depth And recalculate using equation (7) until... The coordinates of the ground projection point can be obtained at that time. Simultaneously, its depth D is obtained; K1 is the camera intrinsic parameter of the tilted image.

[0136] Using the same iterative process, the vertex coordinates of the tilted image are projected onto the ground, resulting in ground point coordinates as follows:

[0137] Then, set the roll and pitch angles of the corrected image to 0, and keep the yaw angle the same as that of the tilt image, i.e., (α,β,γ)=(0,0,γ). corr According to equation (5), the rotation matrix of the corrected image in the world coordinate system is obtained as R2; let the world coordinates of the corrected image in the scene be o. C (X W ,Y W To maintain consistent ground resolution in the central region of the image before and after projection, let h = D. Therefore, the translation vector t2 of the corrected image can be obtained by the following formula.

[0138]

[0139] Finally, use equation (9) to combine the four ground points. By projecting the image forward onto the corrected image, we can obtain the corresponding four image projection points.

[0140] p′=K1[R2p g +t2] (9)

[0141] Substituting the four sets of image coordinates before and after projection into equation (4) yields eight equations. Solving these equations provides the image projection transformation matrix M, which is then applied to equation (3) to transform the tilted image. The transformed spatial coordinates P′(X,Y,Z) are then normalized to image coordinates, i.e. This completes the projection correction of the tilted image.

[0142] The SIFT (Scale Invariant Feature Transform) algorithm is used to extract feature points from the projection-corrected image. Then, the coordinates of the feature points are inversely transformed back to the tilted image (i.e., the original image) using equation (3) to determine the feature points in the tilted image, which are the same as the feature points in the original image. This completes the projection transformation and feature point extraction of the tilted image.

[0143] S3, Matching Image Pair Search:

[0144] Based on feature points in the original image, feature points detected on different images of the same spatial point in a scene are called corresponding points. After extracting these corresponding points, a correspondence needs to be established to estimate the relative pose of the images; this process is called feature point matching. Since large-scale 3D reconstruction tasks involve a large number of images, matching all images is extremely inefficient. Furthermore, because each image covers a small area of ​​the ground, only a few images overlap. Therefore, this invention also utilizes the initial position and pose data of the camera to design a matching image pair search method, significantly reducing the number of image matching operations.

[0145] The basic principle of matching image pairs search is as follows: Figure 5 As shown in the figure, the large dots (actually gray dots) represent the positions of the exposure times of the dual-angle multispectral camera system, and the small dots represent the ground projection points of the center points of the images captured by camera 1 and camera 2, respectively.

[0146] First, the geocentric fixed rectangular coordinates of each original image are obtained using the pose data of each image and equation (6). Simultaneously determine the minimum longitude (lon) in the original image within this dataset. ref latitude ref Let the reference height be alt. ref Taking the spatial point = 0 as the reference point, the coordinates of the reference point for this scene are obtained using equation (6). The coordinates of each original image in the world coordinate system are obtained by equation (10), and the rotation matrix R and translation vector t are solved by equations (5) and (8) respectively.

[0147]

[0148] Then, using the method in equation (7), the center point of each image is projected backwards onto the ground to obtain the projection point O. G (X G ,Y G ,0).

[0149] Let camera 1 be at O C The ground projection coordinates of the center point of a rear-view oblique image taken at the location are O. G Then the neighbor matching image projection point O′ of the reference imageG Satisfying equation (11):

[0150] ||O G -O′ G ||≤r s (11)

[0151] Where, r s The search radius is estimated based on the image overlap rate. Generally, it is sufficient to ensure that the minimum heading and flight strip overlap rates are not less than 50%.

[0152] This search method directly searches for neighboring images with an overlap rate of more than 50% for the reference image for subsequent matching, and the matching efficiency is significantly better than the exhaustive search matching strategy.

[0153] S4. Local Image Rematching:

[0154] To further improve feature point matching accuracy, the homography relationship between image pairs is estimated using the initial matching results, and the feature points of the image to be matched are projected onto the reference image. Based on this, the reference image is divided into blocks, and local rematching is performed within each block. Due to the influence of perspective projection geometric distortion, illumination variation radiation distortion, and other factors, the preceding processing can improve these problems to some extent, but cannot completely solve them. Therefore, the extracted image feature points are not ideal. Thus, local image rematching is used to further improve the effect; however, this step can be omitted.

[0155] In the matching process, two common scenarios, namely mismatch and missed match, are shown in Figures 6(a) and 6(b), leading to poor final matching results. In the high-dimensional space represented by the feature point descriptor, when a non-matching point is the nearest neighbor of the reference point and the matching point is too far from the reference point, a mismatch occurs as shown in Figure 6(a). When there are only a few mismatched pairs between two images, they can be removed using a mismatch elimination algorithm. However, when there are a large number of mismatches, it may lead to errors in subsequent camera pose estimation, thus affecting the accuracy of 3D reconstruction. When the second nearest neighbor of the reference point is too close to the nearest neighbor matching point, the reference point will not match any point because it cannot effectively distinguish between the two points, resulting in a missed match as shown in Figure 6(b). Both of these errors prevent points that could have been correctly matched from being matched, reducing the number of matching point pairs and thus decreasing the robustness of camera pose estimation. To improve these problems, this invention also proposes a local rematching method based on homography transformation guidance. For a pair of images to be matched obtained through back projection search, the feature point matching process is as follows: Figure 7 As shown.

[0156] First, the KD-Tree algorithm from the Fast Library for Approximate Nearest Neighbors (FLANN) is used to search for the nearest and second nearest neighbors of the feature points in the reference image on the matching image. The search is based on the Euclidean distance between the feature point descriptors, as shown in Equation (12), d i d j These are feature point descriptor vectors from the two images, d and d. Euclid It is the Euclidean distance of the vector.

[0157] d Euclid =||d i -d j ||2 (12)

[0158] Let the Euclidean distances between the nearest neighbor, the second nearest neighbor, and the reference point be d1 and d2, respectively. When the two satisfy equation (13), the reference point matches the nearest neighbor. At the same time, r is used as the confidence level of the matched point pair. The smaller the value of r, the higher the accuracy of the matched point pair.

[0159]

[0160] Where η is the matching point filtering threshold, which is generally adjusted according to the required number of matching pairs. The value range is 0-1. The smaller the value, the higher the matching accuracy but the fewer the number of matches.

[0161] For the initial matching process described above, there are generally a large number of mismatched point pairs in the initial matching results. The MAGSAC++ algorithm (an improved edge sample consensus algorithm) is used to remove mismatched point pairs. It selects data points through adaptive geometric sampling and estimates the homography transformation matrix by maximizing the consistency principle, which can achieve a more robust estimation in data with a large amount of noise and outliers.

[0162] The homography transformation matrix M between the two images to be matched is estimated. d Then, referring to equations (3) and (4), the feature point coordinates of the matching image are projected onto the reference image coordinate system (the feature points of the image are the pixels of the image, represented by pixel coordinates; the process of tilting image correction is also a transformation of pixel coordinates, and the two are essentially the same; the implementation process is based on the homography transformation matrix described above and referring to equations (3) and (4)). The reference image is then divided into blocks, with the block method being to divide the reference image into 5×5 sub-images of equal length and width according to the dashed lines, such as... Figure 8 As shown. Then, iterate through each sub-image of the reference image, perform KD-Tree nearest neighbor search matching on the reference feature points and matching feature points within the sub-image region, and then merge the matching results of all sub-images;

[0163] To reduce the impact of projection errors, the reference image and the matching image are swapped (i.e., the reference image is used as the new matching image, and the matching image is used as the new reference image), while simultaneously utilizing the inverse transformation matrix M. d -1 Project the feature points of the new matching image, divide the new reference image into blocks and perform local rematching, and take the intersection of the bidirectional matching results as the matching result of the two images.

[0164] Finally, the MAGSAC++ algorithm is applied to the matching results to remove erroneous matching pairs, resulting in the final matching pairs.

[0165] S5, Multispectral Point Cloud Generation:

[0166] Structure from Motion (SfM) is an algorithm that estimates camera parameters and pose information from multiple overlapping images and reconstructs the 3D structure of a scene. Its core idea is to use images from different viewpoints for feature point matching, estimate and optimize camera parameters and pose based on the spatial geometric constraints of multiple matching image pairs, and then estimate the spatial 3D coordinates of the matching feature points to achieve the reconstruction of a sparse point cloud of a scene. Figure 9 As shown.

[0167] Based on the matching point pairs determined by S4, the precise extrinsic parameters and camera intrinsic parameters (including distortion parameters) of each image are estimated using the open-source 3D reconstruction library OpenSfM. Multi-view dense matching is then performed on the images to obtain depth map point clouds, which form the basis for multispectral point cloud reconstruction. Therefore, using the aforementioned feature point matching results, the image extrinsic parameters and camera intrinsic parameters are estimated through the structure-of-motion reconstruction principle, and then the depth map point cloud of each image is estimated using a multi-view matching method.

[0168] After obtaining the depth map point cloud for each image, the multispectral reflectance map is distorted by combining the camera intrinsic parameters and image extrinsic parameters, and the 5-band spectral values ​​of each point in the depth map point cloud are obtained from it, thus obtaining the multispectral depth map point cloud.

[0169] Then, multispectral depth map point clouds from images acquired by the same camera are stitched together to obtain multispectral depth map point clouds of the same band. In this embodiment, there are 5 sets of point clouds with different spectral bands. Then, the point clouds from two cameras containing different spectral bands are stitched together to generate a complete 10-band multispectral point cloud based on spatial consistency. Finally, outliers are filtered out. The main workflow is as follows: Figure 10 As shown.

[0170] Specifically,

[0171] Since camera lenses generally have a certain degree of radial and tangential distortion, resulting in geometric distortion of the image, in order to improve the accuracy of depth map point cloud spectral reconstruction, it is first necessary to perform distortion correction on the reflectance image. As shown in Equation (14), since the inverse transformation of the distortion model is difficult to solve, the coordinates (x, y) in the distortion-free target image are projected onto the distorted original image through distortion transformation, and the pixel value of the projection position (x', y') is taken as the pixel value of that point in the distortion-free image.

[0172]

[0173] Where r represents the distance from the image point (x,y) to the center of the image plane of the distorted image, and k1, k2, k3 and p1, p2 are the radial and tangential distortion parameters, respectively.

[0174] Based on the image extrinsic parameters and the intrinsic parameters of the two cameras, the forward projection of each spatial point in the depth map point cloud is used to obtain the front pixel coordinates of each spatial point, and then the 5-band spectral reflectance values ​​are obtained from the distortion-corrected reflectance map.

[0175] After obtaining all the multispectral depth map point clouds, these point clouds are divided into two groups, R and B, according to the different cameras used, representing the two multispectral cameras used in this invention. The multispectral depth map point clouds corresponding to the images captured by each camera are stitched and fused to obtain two 5-band multispectral point clouds PC with different spectral bands. R and PC B .

[0176] Then, based on the spatial consistency of corresponding points, a complete 10-band multispectral point cloud is generated from the two 5-band point clouds, such as... Figure 11 As shown, for point cloud PC R Each point in PC B The algorithm searches for its nearest neighbor points. Points with a distance of no more than a threshold d are considered to be of the same name, generating a new spatial point cloud containing the 5-band spectral values ​​of each of the two points. These new spatial points form a 10-band multispectral point cloud. The threshold d is the maximum ground sampling distance of the image, which is adjusted according to the actual situation of the collected data. In this invention, the maximum ground sampling distance of the tilted image is approximately 0.08 meters. Due to the large amount of point cloud data, this invention uses the large-scale vector similarity calculation library hnswlib for nearest neighbor search to improve the algorithm's running speed.

[0177] Finally, to remove point cloud noise, filter out points with spatial anomalies, and make the final point cloud geometry smoother, this invention uses a statistical method from the open-source point cloud data processing library PDAL for anomaly identification. For each point in the multispectral point cloud, k nearest neighbors are searched and the average distance is calculated. And the standard deviation σ, to obtain the outlier judgment threshold shown in equation (15),

[0178]

[0179] Where m is the multiplication factor, and k and m are adjusted appropriately according to the point cloud quality;

[0180] When the distance u of the nearest neighbor of i i Points ≥ t will be marked as outliers and filtered out.

[0181] The three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images proposed in this invention can reconstruct three-dimensional point cloud data from oblique multispectral images containing complex terrain features, improving the representation and analysis capabilities of multidimensional remote sensing data and enabling multispectral stereo remote sensing applications. The advantages of this invention are reflected in the following aspects:

[0182] (a) High Number of Matching Image Feature Points: The multispectral image texture enhancement method in this invention can improve the contrast between vegetation and non-vegetation pixels in the image, enhance the grayscale texture of the image, and reduce the phenomenon of multi-band grayscale nonlinear changes. The tilted image projection correction method can correct large tilt angle images to a downward viewing angle, significantly reducing the degree of perspective geometric distortion. The matching image pair search and local image rematching method can significantly reduce mismatches and increase the number of correct matches. The above processing algorithms greatly improve the feature point extraction and matching logarithm of tilted multispectral images, ensuring the accuracy and integrity of subsequent 3D reconstruction algorithms.

[0183] (b) High integrity of the reconstructed point cloud spatial structure: This invention uses a tilted dual-angle view to acquire multispectral images, which can acquire more object side elevation images compared to a single angle downward view. Combined with projection correction and local rematching methods, more object side elevation reconstructed point clouds can be obtained, and the spatial structure of the reconstructed scene has high integrity.

[0184] (c) Complete Spectral Band Reconstruction of Point Clouds: Since images captured by cameras from different viewpoints cannot be directly registered, existing 3D reconstruction methods cannot directly reconstruct complete spectral band point clouds from unregistered cameras. The depth map-based multispectral point cloud generation method designed in this invention first estimates the depth map point cloud of each image and merges the point clouds of the same camera. Then, it extracts the nearest neighbor points from all camera point clouds to obtain the complete spectral band reflectance value of each spatial point, thereby achieving the generation of complete spectral band multispectral point clouds. Specific Implementation Method Two:

[0186] This embodiment is a three-dimensional multispectral point cloud reconstruction system based on dual-angle multispectral images, including an image texture enhancement unit, a projection correction and feature point extraction unit, a matching image pair search unit, an image matching unit, and a full-band multispectral point cloud generation unit; wherein...

[0187] Image texture enhancement unit: performs texture enhancement on images acquired by different multispectral cameras; the different multispectral cameras are cameras with different spectral bands, and the images acquired by the multispectral cameras are multispectral reflectance maps;

[0188] Projection correction and feature point extraction unit: Based on the texture-enhanced image of each multispectral image, for the tilted image, tilted image projection correction is performed, the SIFT algorithm is used to extract feature points from the projection-corrected image, and then the coordinates of the feature points are inversely transformed back to the tilted image, thereby determining the feature points in the tilted image, which are also the feature points in the original image;

[0189] Matching Image Pair Search Unit: Based on feature points in the original image, the original image is searched for matching image pairs. Neighboring images with an overlap rate exceeding the overlap rate threshold with the reference image are searched for and used for subsequent matching.

[0190] Image matching unit: performs image matching on the searched images;

[0191] Full-band multispectral point cloud generation unit: First, based on the determined matching point pairs, estimate the extrinsic parameters of each image and the intrinsic parameters of the camera, and perform multi-view dense matching on the images to obtain depth map point clouds;

[0192] After obtaining the depth map point cloud for each image, the multispectral reflectance map is distorted by combining the camera intrinsic parameters and image extrinsic parameters, and the multiband spectral values ​​of each point in the depth map point cloud are obtained from it, thus obtaining the multispectral depth map point cloud.

[0193] Then, the multispectral depth map point clouds collected by the same camera are stitched together to obtain multispectral depth map point clouds of the same band. Then, the point clouds containing different bands of spectrum from different multispectral cameras are stitched together to generate a complete full-band multispectral point cloud based on spatial location consistency. Finally, outliers are filtered out.

[0194] The aforementioned full-band multispectral point cloud is a multispectral point cloud covering all bands of different multispectral cameras.

[0195] The system described in this embodiment is a program system corresponding to the three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images. Each processing unit of the system is used to execute the method described in each processing step of the three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images. Specific implementation method three:

[0197] This embodiment is a computer storage medium that stores at least one instruction. The at least one instruction is loaded and executed by a processor to implement the three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images.

[0198] It should be understood that the instructions include computer program products, software, or computerized methods corresponding to any method described in this invention; the instructions can be used to program computer systems or other electronic devices. Computer storage media may include readable media on which instructions are stored, and may include, but are not limited to, magnetic storage media, optical storage media; magneto-optical storage media include read-only memory (ROM), random access memory (RAM), erasable programmable memory (e.g., EPROM and EEPROM), and flash memory layers, or other types of media suitable for storing electronic instructions. Specific implementation method four:

[0200] This embodiment is a three-dimensional multispectral point cloud reconstruction device based on dual-angle multispectral images. The device includes a processor and a memory. It should be understood that this includes any device including a processor and a memory described in this invention. The device may also include other units and modules that perform display, interaction, processing, control, and other functions through signals or instructions.

[0201] The memory stores at least one instruction, which is loaded and executed by the processor to implement the three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images.

[0202] Example

[0203] The three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images in this embodiment is as follows:

[0204] With a flight altitude of 80 meters, a heading, and a flight strip overlap of 80%, dual-angle multispectral images were acquired for a scene with significant differences in the height of ground features and diverse ground feature types. Some images are shown below. Figure 12 As shown, the rear-view tilt camera (hereinafter referred to as "camera 1") acquires 200 5-band images, and the front-view tilt camera (hereinafter referred to as "camera 2") acquires 185 5-band images. The tilt angle of both cameras is 30°. The cameras simultaneously acquire the latitude and longitude coordinates and attitude angle information (i.e., pose data) of the cameras at the exposure time and save them in the image metadata area.

[0205] Each image is processed using the multispectral image texture enhancement method of this invention, and the comparison results are as follows: Figure 13 As shown, Figure 13Figures (a, b, c, and f) show the texture enhancement results and local magnification comparisons. (a) is the original DN value image; (b) is a magnified view of region 1 in (a); (c) is a magnified view of region 2 in (a); (d) is the NDVI-RGB composite image; (e) is a magnified view of region 3 in (d); and (f) is a magnified view of region 4 in (d). A comparison of (a) and (b) shows that the enhanced composite image has higher overall contrast, significantly more geometric features such as points and lines, and a significantly improved texture difference between vegetation and non-vegetation pixels. A comparison of (b) and (e) and (c) and (f) shows that because the spectral reflectance intensity of some man-made features and bare land in this band is similar to that of the surrounding vegetation, these features are difficult to distinguish in the original DN value image. However, in the NDVI-RGB composite image, the multispectral characteristics of the features are utilized, increasing the differences between pixels of different types of features, thus enhancing the image texture. This facilitates the extraction of more and higher-quality feature points, thereby improving the success rate of image matching.

[0206] Applying the tilted image projection correction method of this invention, the comparison results are as follows: Figure 14 As shown, Figure 14 The images show a comparison before and after projection correction, where (a) is the tilted image before correction and (b) is the image after correction. It can be observed that the edges of the building's top surface are not parallel (dashed lines) in the tilted image before correction, but become parallel (dashed lines) after correction. This indicates that projection correction reduces the degree of perspective geometric distortion in the image.

[0207] The matching image pair search method of this invention is used to search for potential matching images of a reference image, and the search results are visualized, such as... Figure 15 As shown, the reference image is IMG_0986_4.tif.

[0208] The SIFT feature point detector is used to extract feature points from the projection-corrected image, and the local image rematching method of this invention is applied to perform feature point matching. The visualization diagram of the matching result is shown below. Figure 16 As shown.

[0209] Table 2 shows a comparison of the feature extraction and matching results with existing open-source frameworks. The results demonstrate that the texture enhancement, projection correction, image pair search, and local re-matching method designed in this invention can obtain a greater number of feature points and significantly increase the number of matching point pairs in the image. Furthermore, it can reduce the reprojection error of matching point pairs to a certain extent, achieving optimal results in all metrics, thus proving the effectiveness of the method designed in this invention.

[0210] Table 2. Image Feature Point Extraction and Matching Results

[0211]

[0212] Applying the multispectral point cloud generation method of this invention to image matching point pairs yields, as shown in the figure. Figure 17 The multispectral dense point cloud model shown is Figure 17 The figures show the visualization results of the multispectral dense point cloud and the vegetation index calculation results. Figure (a) shows the display effect in the visible light true color band, Figure (b) shows the display effect in the standard pseudo-color band, Figure (c) shows the normalized difference vegetation index result of the multispectral point cloud, and Figure (d) shows the green chlorophyll vegetation index result of the multispectral point cloud. The results show that the reconstructed multispectral point cloud spatial structure is relatively complete, and the spectral band recovery is good.

[0213] Applications of extracting individual ground features from point clouds for analyzing their spatial characteristics, such as extracting... Figure 17 (c) Trees and buildings within the Chinese frame, such as Figure 18 As shown. Figure 18 The images show the point cloud extraction results for a single ground feature. (a) shows the normalized difference vegetation index (NDVI) of the point cloud for a single tree. The tree's crown length and width are measured to be 8.27 meters and 7.27 meters, respectively, and its height is 5.76 meters. Combining this with multiple vegetation indices allows for further estimation of biomass, growth stage, and health status. (b) shows the elevation coloring result of the point cloud for a building. Based on the point cloud, the building's roof area is estimated to be approximately 276 square meters, and its volume approximately 2794 cubic meters. The spatial structure information of these ground features is difficult to obtain solely from remote sensing images, and this information is of great significance for applications such as precision agriculture and natural resource monitoring.

[0214] In summary, the three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images designed in this invention can effectively process tilted multispectral images and generate dense point cloud models with relatively complete geometry and spectrum. It overcomes the shortcomings of traditional image-based analysis methods, provides high-quality data for many quantitative remote sensing analyses, and has broad application prospects.

[0215] 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 method for reconstructing three-dimensional multispectral point clouds based on dual-angle multispectral images, characterized in that, Includes the following steps: S1. Perform texture enhancement on images acquired by cameras with different set angles; the camera is a multispectral camera, and the spectral bands corresponding to cameras with different set angles are different, and the images acquired by the multispectral camera are multispectral reflectance maps; Texture enhancement for images acquired by different multispectral cameras includes the following steps: First, images are acquired using different multispectral cameras. Then, normalized difference vegetation index maps are calculated using each camera. ; Define operation rules This indicates that for each element in matrix B, if its value... If the element at the corresponding position in matrix A with the same dimension as B is 1, then the element at the corresponding position is 1; otherwise, it is 0. Represents a numerical value. express , 、<、 One type of symbol; The following formula can be used to remove the shaded areas from the vegetation index map. （1） Vegetation index maps used to represent different multispectral cameras The center wavelength used to represent different multispectral cameras is Single-band reflectivity map; Then, the vegetation pixels in the vegetation index map are combined with the non-vegetation pixels in the blue, green, and red band reflectance maps. Pixels with a vegetation index value less than 0.3 are considered non-vegetation pixels. The calculation formula is as follows: （2） In the formula, These represent the reflectance images for the red, green, and blue bands respectively in the multispectral reflectance diagram; This is the final NDVI-RGB composite image with strong geometric texture; S2. Based on the texture-enhanced image of each multispectral reflectance map, tilt image projection correction is performed for the tilted image. The SIFT algorithm is used to extract feature points from the projection-corrected image, and then the coordinates of the feature points are inversely transformed back into the tilted image to determine the feature points in the tilted image. S3. Based on the feature points in the texture-enhanced image, perform a matching image pair search on the texture-enhanced image, and search for images whose overlap rate with the reference image exceeds the overlap rate threshold for subsequent matching; S4. Perform image matching on the images found in S3; S5. Based on the matching point pairs determined in S4, estimate the extrinsic parameters of each image and the intrinsic parameters of the camera, and perform multi-view dense matching on the images to obtain the depth map point cloud. After obtaining the depth map point cloud for each image, the multispectral reflectance map is distorted by combining the camera intrinsic parameters and image extrinsic parameters, and the multiband spectral values ​​of each point in the depth map point cloud are obtained from it, thus obtaining the multispectral depth map point cloud. Then, the multispectral depth point clouds collected by the same camera are stitched together to obtain multispectral depth point clouds of the same band. Then, based on the consistency of spatial position, the point clouds of different band spectra corresponding to cameras with different settings angles are stitched together to generate a complete full-band multispectral point cloud. Finally, outliers are filtered out. The aforementioned full-band multispectral point cloud is a multispectral point cloud covering all bands of different cameras.

2. The three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images according to claim 1, characterized in that, The process of performing tilted image projection correction includes the following steps: First, determine the perspective projection transformation relationship corresponding to the downward-viewing image after correcting the tilted image. , These represent the image coordinates before perspective projection and the spatial coordinates after projection, respectively. （3） in, It is the image perspective projection matrix. to For elements in the image perspective projection matrix; set up These are the coordinates of the projected image. Equation (3) can be used to establish the equation shown in equation (4). （4） make Then the number of unknown variables in the equation becomes 8, and the matrix can be solved by finding 4 sets of image coordinates before and after projection. Use the coordinates of the four vertices of the tilted image as... Solve for the corresponding 4 projection points : First, determine the homogeneous pixel coordinates of the center point of the tilted image. Projected backwards onto a horizontal ground, let the world coordinates of the projection point be... The rotation matrix of the tilted image in the world coordinate system can be obtained from the pose data. ; Then, the latitude, longitude, and altitude coordinates of each camera exposure position are converted to geocentric fixed rectangular coordinates to obtain the translation vector of the image in the world coordinate system. ; Based on the rotation matrix of the image in the world coordinate system The ground projection point of the image center point is obtained according to equation (7): （7） In the formula, For the camera's intrinsic parameters of the tilted image; This represents the height from the camera's optical center to the ground, i.e., the depth of the projection point; let the initial depth be... For the drone's flight altitude, when Time-updated depth And recalculate using equation (7), The world coordinate Z-axis reference threshold, until The coordinates of the ground projection point can be obtained at that time. At the same time, obtain depth ; Projecting the vertex coordinates of the tilted image onto the ground yields the coordinates of the ground points. ; Then, set the roll and pitch angles of the corrected image to 0, and keep the yaw angle the same as that of the tilt image. And the rotation matrix of the corrected image in the world coordinate system is obtained as follows. ; Let the corrected image be in world coordinates as ,make The translation vector of the corrected image is obtained. ; （8） Finally, use equation (9) to combine the four ground points. By projecting the image forward onto the corrected image, we can obtain the corresponding four image projection points. （9） Substituting the image coordinates before and after the four projections into equation (4) yields eight equations. Solving these equations provides the matrix. Then, equation (3) is applied to transform the tilted image, and the transformed spatial coordinates are... Normalized to image coordinates, i.e. This completes the projection correction of the tilted image.

3. The method for reconstructing three-dimensional multispectral point clouds based on dual-angle multispectral images according to claim 2, characterized in that, The process of performing image pair matching search on a texture-enhanced image based on feature points in the texture-enhanced image includes the following steps: First, the geocentric fixed Cartesian coordinates of each texture-enhanced image are obtained based on its pose data. Simultaneously determine the minimum longitude in the texture-enhanced image dataset. ,latitude Let the reference height be Using a spatial point as a reference point, the geocentric fixed rectangular coordinates of the reference point are obtained as follows: The coordinates of each texture-enhanced image in the world coordinate system are obtained from equation (10), and the rotation matrix in the world coordinate system is solved separately. Translation vector ; Then, the center point of each image is projected backwards onto the ground to obtain the projection point. ; Suppose a certain camera is in The ground projection coordinates of the center point of a rear-view oblique image taken at the location are: Then the neighboring image projection points of the reference image are matched. Satisfying equation (11): ||O G -O′ G ||≤r s (11) in, The search radius is estimated based on the image overlap rate. Images whose overlap with the reference image exceeds an overlap threshold are searched for for subsequent matching.

4. The three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images according to claim 3, characterized in that, The process of image matching for images found by S3 includes the following steps: First, the KD-Tree algorithm is used to search for the nearest and second nearest neighbors of feature points in the reference image on the matching image. Let the Euclidean distances between the nearest and second nearest neighbors and the feature points in the reference image be respectively... , When both satisfy equation (13), the feature points in the reference image are matched with their nearest neighbors, and at the same time, the feature points in the reference image are matched with their nearest neighbors. As the confidence level of the matched point pair; （13） in, It is the matching point filtering threshold; For the matching results, the MAGSAC++ algorithm is used to remove mismatched point pairs. In the process of removing mismatched point pairs, the homography transformation matrix is ​​estimated by maximizing the consistency principle. The homography transformation matrix between the two images to be matched is obtained. Next, the coordinates of the feature points in the matching image are projected onto the coordinate system of the reference image, and the reference image is divided into blocks, each block being a sub-image. Then, each sub-image of the reference image is traversed, and local re-matching is performed, that is, KD-Tree nearest neighbor search matching is performed on the reference feature points and matching feature points within the sub-image region, and then the matching results of all sub-images are merged. Swap the reference image and the matching image, and simultaneously utilize Project the feature points of the new matching image, divide the new reference image into blocks, and perform local rematching. The intersection of the bidirectional matching results is taken as the matching result of the two images; Finally, the MAGSAC++ algorithm is used to remove mismatched point pairs from the matching results to obtain the final matching point pairs.

5. The three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images according to claim 4, characterized in that, In S5, the process of estimating the extrinsic parameters of each image and the intrinsic parameters of the camera based on the matching point pairs determined in S4 is achieved through the principle of motion reconstruction structure.

6. The method for reconstructing three-dimensional multispectral point clouds based on dual-angle multispectral images according to claim 5, characterized in that, In the final outlier filtering process described in S5, outlier identification is performed using a statistical method from the open-source point cloud data processing library PDAL. Search for each point in the multispectral point cloud Find the nearest neighbors and calculate the average distance. and standard deviation The outlier detection threshold is obtained as follows: （15） in, It is the multiplication factor; When the Distance between nearest neighbors Points that are not found will be marked as outliers and filtered out.

7. A three-dimensional multispectral point cloud reconstruction system based on dual-angle multispectral images, characterized in that, It includes an image texture enhancement unit, a projection correction and feature point extraction unit, a matching image pair search unit, an image matching unit, and a full-band multispectral point cloud generation unit; among which, Image texture enhancement unit: performs texture enhancement on images acquired by different multispectral cameras; the different multispectral cameras are cameras with different spectral bands, and the images acquired by the multispectral cameras are multispectral reflectance maps; Texture enhancement for images acquired by different multispectral cameras includes the following steps: First, images are acquired using different multispectral cameras. Then, normalized difference vegetation index maps are calculated using each camera. ; Define operation rules This indicates that for each element in matrix B, if its value... If the element at the corresponding position in matrix A with the same dimension as B is 1, then the element at the corresponding position is 1; otherwise, it is 0. Represents a numerical value. express , 、<、 One type of symbol; The following formula can be used to remove the shaded areas from the vegetation index map. （1） Vegetation index maps used to represent different multispectral cameras The center wavelength used to represent different multispectral cameras is Single-band reflectivity map; Then, the vegetation pixels in the vegetation index map are combined with the non-vegetation pixels in the blue, green, and red band reflectance maps. Pixels with a vegetation index value less than 0.3 are considered non-vegetation pixels. The calculation formula is as follows: （2） In the formula, These represent the reflectance images for the red, green, and blue bands respectively in the multispectral reflectance diagram; This is the final NDVI-RGB composite image with strong geometric texture; Projection correction and feature point extraction unit: Based on the image after texture enhancement of each multispectral reflectance map, for the tilted image, tilted image projection correction is performed, the SIFT algorithm is used to extract feature points from the projection-corrected image, and then the coordinates of the feature points are inversely transformed back to the tilted image to determine the feature points in the tilted image. Matching Image Pair Search Unit: Based on feature points in the texture-enhanced image, the matching image pair search is performed on the texture-enhanced image, and neighboring images with an overlap rate exceeding the overlap rate threshold with the reference image are searched for subsequent matching; Image matching unit: performs image matching on the searched images; Full-band multispectral point cloud generation unit: First, based on the determined matching point pairs, estimate the extrinsic parameters of each image and the intrinsic parameters of the camera, and perform multi-view dense matching on the images to obtain depth map point clouds; After obtaining the depth map point cloud for each image, the multispectral reflectance map is distorted by combining the camera intrinsic parameters and image extrinsic parameters, and the multiband spectral values ​​of each point in the depth map point cloud are obtained from it, thus obtaining the multispectral depth map point cloud. Then, the multispectral depth point clouds collected by the same camera are stitched together to obtain multispectral depth point clouds of the same band. Then, based on the consistency of spatial position, the point clouds of different band spectra corresponding to cameras with different settings angles are stitched together to generate a complete full-band multispectral point cloud. Finally, outliers are filtered out. The aforementioned full-band multispectral point cloud is a multispectral point cloud covering all bands of different multispectral cameras.

8. A computer storage medium, characterized in that, The storage medium stores at least one instruction, which is loaded and executed by a processor to implement the three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images as described in any one of claims 1 to 6.

9. A three-dimensional multispectral point cloud reconstruction device based on dual-angle multispectral images, characterized in that, The device includes a processor and a memory, the memory storing at least one instruction, which is loaded and executed by the processor to implement the three-dimensional multispectral point cloud reconstruction method based on dual-angle multispectral images as described in any one of claims 1 to 6.

Images