A method for stitching UAV hyperspectral images based on shallow water bottom feature matching

Abstract

The invention discloses a method for stitching unmanned aerial vehicle hyperspectral images based on shallow water bottom feature matching, which belongs to the field of image data processing technology, and includes: obtaining unmanned aerial vehicle hyperspectral images to be stitched; flare data correction step; inverting the bottom information and water depth corresponding to each pixel of each unmanned aerial vehicle hyperspectral image; geometrically correcting the unmanned aerial vehicle hyperspectral image according to the longitude and latitude information; unifying the resolution of the images to be stitched; selecting the correction image sequence; calculating the overlapping area between the registered images; extracting the key features of the land and water areas; obtaining the matching point pair array; constructing a transformation model; calculating the new coordinates of the image, and stitching the two adjacent images to be stitched; and comprehensively fusing the hyperspectral data, the water depth data and the bottom data to output the final stitched image. The invention utilizes the characteristic that the seabed topography remains unchanged, inverts the bottom topography through the hyperspectral, and identifies the features according to the seabed information, thereby improving the stitching accuracy.

Description

Technical Field

[0001] The invention belongs to the technical field of remote sensing image processing, and in particular relates to a method for stitching unmanned aerial vehicle hyperspectral images based on shallow water bottom feature matching. Background Art

[0002] In the context of today's rapid technological development, efficient, flexible and high-resolution data acquisition has become increasingly urgent for ocean research and management. In this field, two different remote sensing platforms, satellite remote sensing and drone remote sensing, have emerged, each with its own unique characteristics and advantages. Space-based remote sensing uses high-altitude platforms such as satellites to achieve global macroscopic observations, while air-based remote sensing uses low-altitude platforms such as drones to achieve high-resolution monitoring of local areas in a more flexible way.

[0003] UAV remote sensing, through the application of low-altitude platforms such as drones, has demonstrated unique advantages in terms of high resolution, immediacy and flexibility. Especially in the marine field, the demand for UAV remote sensing has gradually emerged. It can respond quickly, monitor in real time and obtain data with high precision, which is more suitable for local needs such as marine ecological monitoring, resource management, and border security.

[0004] The single remote sensing image acquired by the UAV remote sensing platform has a small coverage area. For a large area, multiple remote sensing images need to be taken for image stitching and fusion. The stitching technology based on image feature points is a relatively mature image stitching technology currently studied. This technology achieves image stitching by finding points in the image where the regional energy changes faster and matching the feature points of the same name in different images.

[0005] The stitching of marine remote sensing images has always been a difficult problem. It is difficult to set up control points on the sea surface, and feature points are scarce. Therefore, the traditional registration stitching based on control points and image feature matching is not suitable for marine areas. The image correction technology based on the POS system is too expensive. A high-precision POS system often costs tens of thousands or even hundreds of thousands, which limits its application. The traditional image stitching technology based on feature points is generally based on band grayscale recognition, which is very difficult for the ocean, and for shallow water areas, the bottom information is covered.

[0006] Due to the large number of remote sensing images of large-scale ocean waters, the objects collected are relatively simple and uniform, and the use of this technology has problems such as few feature points and large registration errors, resulting in the stitching effect cannot meet the application requirements of marine drone remote sensing, and large-scale image stitching often fails. Image stitching using geographic coordinate information requires high-precision GPS information, otherwise the geographic information obtained from remote sensing images has large errors, resulting in positional misalignment between stitched images, and the stability and wind resistance of drones and the performance of GPS measurement equipment make it difficult to guarantee the accuracy of the POS information corresponding to each image.

[0007] Chinese patent CN114936971A combines the existing remote sensing image stitching technology and the characteristics of multispectral images in fishery waters, and proposes a method and system for stitching multispectral images of drone remote sensing in waters. This method is suitable for remote sensing positioning and image stitching in fishery waters. The research area is small and it is not suitable for large ocean fields. There are a large number of grids in its test conditions that can be used for feature matching, and it does not solve the problem of pure water stitching, and is not suitable for most marine drone remote sensing applications.

[0008] The paper "Research on Processing and Stitching Algorithms of UAV Ocean Remote Sensing Images [D]. Yanshan University, 2022. DOI: 10.27440 / d.cnki.gysdu.2022.001608. " conducted an in-depth study of the image preprocessing and stitching algorithms of ocean remote sensing images, and provided a stitching process and algorithm suitable for UAV ocean remote sensing images. This method is used for UAV ocean image stitching, but it only optimizes the feature matching and projection algorithms, and does not optimize the characteristics of ocean remote sensing images. Its research method is also only applicable to shore, where the water is slightly deeper, and the water body blurs the original characteristic information of the bottom. In some pure water images, the feature points are basically undetectable, resulting in large limitations in application scenarios. Summary of the invention

[0009] In order to solve the technical problems that the existing UAV land remote sensing image stitching method is not suitable for ocean remote sensing image stitching, and the research area is small and not suitable for a large ocean field of view, the present invention proposes a UAV hyperspectral image stitching method based on shallow water bottom feature matching, which can solve the above problems.

[0010] In order to solve the above technical problems, the present invention adopts the following technical solutions to achieve the above problems:

[0011] A method for stitching UAV hyperspectral images based on shallow water bottom feature matching includes the following steps:

[0012] Step 1: Obtain the drone hyperspectral images to be stitched, and simultaneously obtain the geographic coordinate projection information corresponding to each image and the latitude and longitude corresponding to each image;

[0013] Step 2, flare data correction step, calculates the radiation value of each UAV hyperspectral image, and uses the radiation value to correct the flare data in the UAV hyperspectral image;

[0014] Step 3: Invert the bottom information and water depth corresponding to each pixel of each UAV hyperspectral image, and generate a matrix with the same size as the original image;

[0015] Step 4: geometrically correct the UAV hyperspectral image according to the latitude and longitude information to obtain a corrected image, and add geographic coordinate projection information to the corrected image;

[0016] Step 5: resample all the images to be stitched to a uniform resolution;

[0017] Step 6: Select a corrected image sequence, perform rough image registration according to the geographic coordinate projection information of each corrected image, and obtain a registered image of the corrected image and a translation transformation matrix between adjacent registered images;

[0018] Step 7, calculating the overlapping area between the registered images, and extracting the overlapping area image to obtain the image to be stitched;

[0019] Step 8, using deep learning technology to automatically identify and extract key features of land and water areas, separate the land area and water area in the image to be stitched, and process the land area and water area respectively to obtain feature points of each land feature band image, as well as feature points of each water feature band image, feature points of the bottom information matrix, and feature points of the water depth;

[0020] Step nine, respectively determine the total number of feature points of each feature band image of the land area and the total number of feature points of each feature band image of the water area; when the total number is not less than N, first identify and delete the multi-band repeated feature points; avoid repeated matching, establish a hypothetical matching set of feature points between two adjacent images to be spliced, and use spatial distance as a similarity evaluation index of feature points, screen out matching points whose spatial distance meets the set value, and obtain a matching point pair array;

[0021] Step ten, constructing a transformation model according to the matching point pair array;

[0022] Step 11, calculating the new coordinates of the image according to the transformation model, and splicing the two adjacent images to be spliced;

[0023] Step 12: After all the images to be stitched are stitched together, the hyperspectral data, water depth data and bottom data are integrated to output the final stitched image.

[0024] In some embodiments, the step eight includes:

[0025] For the land area, the image of each characteristic band is extracted respectively to obtain the land characteristic band image, and the accelerated robust feature model SURF1 is constructed for each land characteristic band image to obtain the feature points of each land characteristic band image; for the water body area, the image of each characteristic band is extracted respectively to obtain the water body characteristic band image, and the bottom sediment information matrix and the water depth matrix of the water body area are extracted, and the accelerated robust feature model SURF2 is constructed for each water body characteristic band image to obtain the feature points of each water body characteristic band image, the feature points of the bottom sediment information matrix and the feature points of the water depth.

[0026] In some embodiments, for water body areas, the total number of feature points of the bottom information matrix and the total number of feature points of the water depth matrix are also determined. When the total number of feature points of the water body characteristic band image, the total number of feature points of the bottom information matrix and the total number of feature points of the water depth matrix are all less than N, the transformation model cannot be constructed; if the total number of feature points of the water body characteristic band image is greater than N, and the total number of feature points of the bottom information matrix and the total number of feature points of the water depth matrix are less than N, the transformation model is generated according to the matching point pair array; when none of the above situations are met, the feature points of the bottom information matrix and the feature points of the water depth matrix are queried and indexed by a KD tree method to establish a hypothetical matching set of feature points between two adjacent images to be spliced, and the spatial distance is used as a similarity evaluation index of the feature points, and matching points whose spatial distance meets the set value are screened out to obtain a second matching point pair array, and the distance histogram is used to constrain the global similarity attribute to eliminate mismatched point pairs, update the matching point pair array, and generate a transformation model according to the second matching point pair array.

[0027] In some embodiments, step 2 includes:

[0028] The calculation method of the radiation value is:

[0029] L(λ)=F L (λ)×(DN r (λ)-DN amb (λ));

[0030] Where λ is the wavelength, F L (λ) is the calibration coefficient, DN r (λ) is the brightness value of the pixel, DN amb (λ) is noise;

[0031] Methods for correcting flare data in UAV hyperspectral images using radiation values ​​include:

[0032] The flare data is identified using a convolutional neural network and classified according to the strength of the flare. The overexposed flare areas are masked, and the remaining flare areas are divided into strong flare areas and weak flare areas. The strong flare areas and weak flare areas are corrected separately using surrounding pixels.

[0033] In some embodiments, step one further includes obtaining the coordinates of the projection center in the ground photogrammetry coordinate system, and the geometric correction in step four includes:

[0034] (41) Establish an image space auxiliary coordinate system parallel to the ground photogrammetry system;

[0035] (42) The image space auxiliary coordinates (X, Y, Z) of the pixel and the ground photogrammetric coordinates (X A ,Y A , Z A ) between:

[0036] ;

[0037] Among them, X S , Y S , Z S are the coordinates of the pixel projection center in the ground photogrammetry coordinate system, and α is the scale factor.

[0038] In some embodiments, adding geographic coordinate projection information to the corrected image includes:

[0039] The remote sensing image coordinates and geographic coordinates are converted using affine matrix parameters, which have 6 parameters, namely ,in, , represents the geographic projection coordinates of the upper left pixel of the corrected image, , Respectively represent the ground resolution of the corrected image pixels in longitude and latitude, Represents the sine value of the rotation angle of the corrected image;

[0040] Rotate the corrected image clockwise around the center 0 point Degrees, and then calculate the coordinates of point 0 based on the corresponding longitude and latitude. , , The calculation formula is:

[0041] ;

[0042] in and Respectively represent the width and height of the corrected image;

[0043] Correction image , Any pixel other than Geographic projection coordinates of The calculation formula is:

[0044] .

[0045] In some embodiments, the method of rough image registration in step 6 is:

[0046] Two adjacent corrected images and images The geographic coordinates of the upper left corner are and , to correct the image As the benchmark, the corrected image is calculated Relative to the calibrated image The offset:

[0047] ;

[0048] Correct the image and images Separately around and Rotation yaw angle Degrees to the route direction, correct the image The pixels in the corrected image can be calculated by the translation transformation matrix The image coordinates in .

[0049] In some embodiments, in step seven, the method for calculating the overlapping area between the registered images includes:

[0050] Calculate the overlapping area of ​​the quadrilateral formed by the four corners of the two registered images respectively;

[0051] Add a certain offset δ outward from the overlapping area to obtain the coordinates of the overlapping area;

[0052] The coordinates of the overlapping area are used to create a mask image, and an image containing only the overlapping area is obtained from the registered image, which are respectively recorded as image and .

[0053] In some embodiments, in step nine, the method for screening matching points includes:

[0054] Using a query index KD tree method, preliminarily searching for n feature points in the second image to be stitched that match each feature point in the first image to be stitched, forming a one-to-many hypothetical matching set, wherein the first image to be stitched and the second image to be stitched are two images to be stitched;

[0055] Calculate the spatial distance between each feature point in the first image to be stitched and the feature points in the assumed matching set;

[0056] Find the nearest neighbor spatial distance and the next nearest neighbor spatial distance, and calculate the ratio of the nearest neighbor spatial distance to the next nearest neighbor spatial distance ,when When , the two feature points are matched point pairs, and all matching point pairs are found to obtain a matching point pair array.

[0057] In some embodiments, constructing a transformation model according to the matching point pair array in step 10 includes:

[0058] For water areas, when the total number of feature points of the corresponding characteristic band images is not less than N, and the total number of feature points of the bottom information matrix and the total number of feature points of the water depth matrix are both less than N, a transformation model is generated based on the matching point pair array; otherwise, the distance histogram is used to constrain the global similarity attribute to eliminate erroneous matching point pairs, the matching point pair array is updated, and a transformation model is generated based on the updated matching point pair array.

[0059] Compared with the prior art, the advantages and positive effects of the present invention are as follows: the UAV hyperspectral image stitching method based on shallow water bottom feature matching of the present invention utilizes the characteristic of unchanged seabed topography, inverts the seabed topography through hyperspectral, identifies features based on seabed information, and thereby improves stitching accuracy, and is particularly suitable for stitching of shallow water ocean hyperspectral images.

[0060] Other features and advantages of the present invention will become more apparent after reading the detailed description of the embodiments of the present invention in conjunction with the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0062] Figure 1 It is a flow chart of an embodiment of the method for stitching UAV hyperspectral images based on shallow water bottom feature matching proposed by the present invention;

[0063] Figure 2 It is a UAV hyperspectral image in one embodiment of the UAV hyperspectral image stitching method based on shallow water bottom feature matching proposed by the present invention;

[0064] Figure 3 is based on Figure 2 Inverted bathymetric map;

[0065] Figure 4 Yes Figure 2 Geometrically corrected image;

[0066] Figure 5 It is a stitched image in an embodiment of the UAV hyperspectral image stitching method based on shallow water bottom feature matching proposed by the present invention. DETAILED DESCRIPTION

[0067] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0068] Example 1, see Figure 1 As shown, this embodiment proposes a method for stitching UAV hyperspectral images based on shallow water bottom feature matching, including the following steps:

[0069] Step 1: Obtain the drone hyperspectral images to be stitched, and at the same time obtain the geographic coordinate projection information corresponding to each image and the longitude and latitude corresponding to each image.

[0070] When the drone is in flight, the hyperspectral images obtained usually carry supporting POS data. This makes it easier to process images. Receive the hyperspectral images of drone remote sensing. The original data exists in raw format and is decoded and synthesized into .nc format. And according to the longitude and latitude POS information corresponding to the single-band image, the abnormal values ​​of the longitude and latitude data at takeoff and return are eliminated; the longitude and latitude POS information includes longitude lon, latitude lat, altitude H, heading angle κ, pitch angle φ, roll angle ω, ground resolution resolution, and the geographic coordinate system and projection coordinate system of the measurement system; the calculation formula of ground resolution E is as follows; resolution=pixelsize*H / f.

[0071] IMU data mainly includes three data: heading angle (Phi), pitch angle (Omega) and roll angle (Kappa).

[0072] GPS data is generally expressed as X, Y, and Z, representing the geographic location of the aircraft at the time of exposure during flight.

[0073] Step 2, the flare data correction step, calculates the radiation value of each UAV hyperspectral image and uses the radiation value to correct the flare data in the UAV hyperspectral image.

[0074] Adaptive solar flare correction technology is used to dynamically adjust mask parameters according to the light intensity and angle during image acquisition, effectively reducing image distortion caused by solar reflection, classifying flares, correcting weak flares in the image, and correcting or masking strong flares according to their characteristics.

[0075] Step three: Invert the bottom information and water depth corresponding to each pixel of each UAV hyperspectral image and generate a matrix with the same size as the original image.

[0076] Step 4: geometrically correct the UAV hyperspectral image according to the latitude and longitude information to obtain the corrected image, and add geographic coordinate projection information to the corrected image. This step is used to establish the spatial mapping relationship between pixels and corresponding ground points.

[0077] Step 5: resample all the images to be stitched to a uniform resolution.

[0078] Step six, select the correction image sequence, perform rough image registration according to the geographic coordinate projection information of each correction image, and obtain the registration image of the correction image and the translation transformation matrix between adjacent registration images.

[0079] Step seven, calculate the overlapping area between the registered images, and extract the overlapping area image to obtain the image to be stitched.

[0080] Step eight, use deep learning technology to automatically identify and extract key features of land and water areas, separate the land areas and water areas in the images to be stitched, and process the land areas and water areas separately to obtain the feature points of each land feature band image, as well as the feature points of each water feature band image, the feature points of the bottom information matrix, and the feature points of the water depth.

[0081] Step nine, respectively determine the total number of feature points of each feature band image in the land area and the total number of feature points of each feature band image in the water area; when the total number is not less than N, first identify and delete the multi-band repeated feature points; avoid repeated matching, establish a hypothetical matching set of feature points between two adjacent images to be stitched, and use spatial distance as the similarity evaluation index of feature points to screen out matching points whose spatial distance meets the set value to obtain a matching point pair array.

[0082] Step 10: construct a transformation model based on the matching point pair array.

[0083] Step eleven: calculate the new coordinates of the image according to the transformation model, and splice the two adjacent images to be spliced.

[0084] Step 12: After all the images to be stitched are stitched together, the hyperspectral data, water depth data and bottom data are integrated to output the final stitched image.

[0085] like Figure 2 As shown in the figure, the bottom features are not obvious due to the influence of water. Figure 3 As shown, it is an inverted water depth image. It can be seen that the inverted water depth map has obvious features and is suitable for feature matching. The UAV hyperspectral image stitching method based on shallow water bottom feature matching of this embodiment, 1. In view of the fact that there is less detailed information in marine remote sensing images, the existing UAV remote sensing image stitching technology is all for land, and there is a lack of pure water remote sensing image stitching methods. This technology is based on the direct stitching and coarse registration based on the POS system. A method for stitching UAV hyperspectral images based on shallow water characteristics is designed, and the seabed sediment information and water depth information are introduced. In the sea area far from the shore, without the help of control points, it provides an effective way for the precise registration of UAV images, and improves the accuracy of marine UAV remote sensing operations within 30m water depth.

[0086] 2. The traditional method of feature detection based on low-dimensional grayscale images and multi-dimensional reduced-dimensional images has been improved. The advantages of high spectral resolution of hyperspectral sensors have been fully utilized to detect more feature points and achieve higher image stitching accuracy. Some characteristic bands that are sensitive to water bodies are used for matching, which is more effective than traditional RGB bands or one-dimensional grayscale images. By comparing and finding bands suitable for matching feature points of offshore shallow water images, the matching efficiency is improved.

[0087] Considering the differences in key features between land and water areas, in some embodiments, step eight further includes:

[0088] For land areas, images of each characteristic band are extracted respectively to obtain land characteristic band images, and accelerated robust feature models SURF1 are constructed for each land characteristic band image to obtain the feature points of each land characteristic band image; for water areas, images of each characteristic band are extracted respectively to obtain water characteristic band images, and the bottom sediment information matrix and water depth matrix of the water area are extracted, and accelerated robust feature models SURF2 are constructed for each water characteristic band image to obtain the feature points of each water characteristic band image, the feature points of the bottom sediment information matrix and the feature points of the water depth.

[0089] In some embodiments, deep learning technology can be used to automatically identify and extract key features of land and water areas. By extracting key features of land and water areas respectively, the extracted feature points are more accurate, thereby improving recognition accuracy.

[0090] In some embodiments, for water areas, it also includes determining the total number of feature points of the bottom information matrix and the total number of feature points of the water depth matrix. When the total number of feature points of the water body characteristic band image, the total number of feature points of the bottom information matrix and the total number of feature points of the water depth matrix are all less than N, the transformation model cannot be constructed; if the total number of feature points of the water body characteristic band image is greater than N, and the total number of feature points of the bottom information matrix and the total number of feature points of the water depth matrix are less than N, the transformation model is generated according to the matching point pair array; when none of the above situations are met, the query index KD tree method is used for the feature points of the bottom information matrix and the feature points of the water depth matrix to establish a hypothetical matching set of feature points between two adjacent images to be spliced, and the spatial distance is used as the similarity evaluation index of the feature points, and the matching points whose spatial distance meets the set value are screened out to obtain a second matching point pair array, and the distance histogram is used to constrain the global similarity attribute to eliminate mismatched point pairs, update the matching point pair array, and generate a transformation model according to the second matching point pair array.

[0091] In some embodiments, step 2 includes:

[0092] The calculation method of the radiation value L(λ) is:

[0093] L(λ)=F L (λ)×(DN r (λ)-DN amb (λ));

[0094] Where λ is the wavelength, F L (λ) is the calibration coefficient, DN r (λ) is the brightness value of the pixel, DN amb (λ) is noise, which is a known value or can be obtained from the image.

[0095] In step 2, the method of using the radiation value to correct the flare data in the UAV hyperspectral image includes:

[0096] Convolutional neural network recognition is used to identify flare data and classify them according to the strength of the flare. The overexposed flare areas are masked, and the remaining flare areas are divided into strong flare areas and weak flare areas. The surrounding pixels are used to correct the strong flare areas and weak flare areas respectively.

[0097] In this scheme, the characteristic that the infrared band radiation in the water area tends to 0 is used to identify flares in the ocean. When correcting the strong flare areas and the weak flare areas, a statistical method is used to estimate the original signal, and the statistical characteristics of the non-flare area are used to estimate the true reflectivity of the area covered by the flare.

[0098] In some embodiments, in step three, the method for inverting water depth and bottom information is as follows:

[0099] Based on the radiation transfer equation, the shallow water hyperspectral remote sensing reflectance can be expressed as:

[0100] .

[0101] In the formula, is the remote sensing reflectance just below the surface; and is a constant; and are the sun and sensor zenith angles, respectively; is the vertical attenuation coefficient of water body; is the seafloor reflectivity; is the water depth. In the above formula, except for the water depth Except for this, all other parameters are functions of wavelength.

[0102] Different substrate types (such as sand, mud, rock, coral, etc.) have different normalized substrate indexes (NBIs). The random forest method is used to classify the substrate types, determine the number of trees, the maximum depth of the tree, the minimum number of samples required for each node split, and other parameters. Cross-validation (such as k-fold cross-validation) is used to evaluate the generalization ability of the model and find the best parameter configuration. The random forest model is trained using the selected features and configurations, and the trained model is applied to the new hyperspectral dataset for substrate classification. Water depth inversion is performed according to different substrate types.

[0103] The formula can be expressed as The function of 6 variables is as follows:

[0104] .

[0105] Hyperspectral sensors have far more than 6 spectral channels. The solution process usually uses a nonlinear optimization algorithm, that is, first establish a UAV observation And formula simulation The error function between them is then used to adjust the values ​​of the six variables through a nonlinear optimization algorithm. When the error function reaches a minimum, it is the solution to the equation.

[0106] Solve the water depth and bottom information corresponding to each pixel and generate a matrix of the same size as the original image.

[0107] In some embodiments, step 1 also includes obtaining the coordinates of the projection center in the ground photogrammetry coordinate system. In order to establish a spatial mapping relationship between pixels and corresponding ground points, an image space auxiliary coordinate system parallel to the ground photogrammetry system must be established first. If the coordinates of pixel a in the image space auxiliary coordinate system and the image space coordinate system are (X, Y, Z) and (x, y, -f), respectively, the coordinates of the corresponding ground point A and the projection center S in the ground photogrammetry coordinate system are (X, Y, Z) and (x, y, -f). A , Y A , Z A ) and (X S , Y S , Z S ), since the corresponding axes of the ground photogrammetry coordinate system and the image space auxiliary coordinate system are parallel to each other, according to the similar triangle relationship, the image space auxiliary coordinates (X, Y, Z) of the pixel and the ground photogrammetry coordinates (X A , Y A , Z A ) between them.

[0108] Therefore, the geometric correction in step 4 includes:

[0109] (41) Establish an image space auxiliary coordinate system parallel to the ground photogrammetry system;

[0110] (42) The image space auxiliary coordinates (X, Y, Z) of the pixel and the ground photogrammetric coordinates (X A ,Y A , Z A ) between:

[0111] ;

[0112] Among them, X S , Y S , Z S are the coordinates of the pixel projection center in the ground photogrammetry coordinate system, and α is the scale factor.

[0113] The inverse formula of the collinearity condition equation is:

[0114] .

[0115] Where f is the focal length of the camera, , is the coordinate of the pixel in the image, is the coordinate of the projection center S in the ground photogrammetry coordinate system, which can be obtained from the POS system (or GPS system). , , is the direction cosine, which is the cosine of the angle between each axis of the image space coordinate system and each axis of the corresponding image space auxiliary coordinate system. , Roll Angle , Yaw angle The sine and cosine components are shown in the following formula:

[0116] .

[0117] That is, the spatial coordinates of the ground point corresponding to the pixel can be solved.

[0118] According to the above correction mathematical model, the new coordinates of the four corners of the image after correction are calculated, the transformation matrix is ​​obtained through the old coordinates and the new coordinates, and the new corrected image is calculated using the bilinear interpolation resampling method to complete the geometric correction of the image and update the geographic coordinate information of the image.

[0119] like Figure 2 , Figure 4 As shown in the figure, it is a comparison diagram before and after geometric correction. Figure 4 is the corrected image.

[0120] In some embodiments, in step 4, adding geographic coordinate projection information to the corrected image includes:

[0121] The affine matrix parameters are used to transform the remote sensing image coordinates and geographic coordinates. There are 6 affine matrix parameters, which are: ,in, , represents the geographic projection coordinates of the upper left pixel of the corrected image, , Respectively represent the ground resolution of the corrected image pixels in longitude and latitude, Represents the sine of the rotation angle of the corrected image.

[0122] Rotate the corrected image clockwise around the center 0 point Degrees, and then calculate the coordinates of point 0 based on the corresponding longitude and latitude. , , The calculation formula is:

[0123] ;

[0124] in and Represent the width and height of the corrected image respectively.

[0125] Correction image , Any pixel other than Geographic projection coordinates of The calculation formula is:

[0126] .

[0127] In some embodiments, the resampling method uses bilinear interpolation. If the row and column coordinates of the corresponding image point P on the corrected image on the original image are (i, j), then take the image point to the right of the image point (i + 1, j), the image point above (i, j + 1), and the image point at the upper right corner (i + 1, j + 1). The grayscale values ​​corresponding to the above four points are represented by f(i, j), f(i +1, j), f(i, j + 1), and f(i + 1, j + 1), respectively, and are represented by dx and dy, that is, the decimal part of the corresponding coordinates. Finally, a weighted operation is performed according to the distance of the subscript value from the surrounding points to obtain the grayscale value fp that should be assigned to the orthophoto. The specific calculation formula is as follows:

[0128] fp=(1-dx)(1-dy)f(i,j)+(1-dx)dyf(i+1,j)+dx(1-dy)f(i,j+1)+dxdyf(i+1,j+1).

[0129] In some embodiments, the method of rough image registration in step 6 is:

[0130] Assume that two adjacent corrected images on the same route and images The geographic coordinates of the upper left corner are and , to correct the image As the benchmark, the corrected image is calculated Relative to the calibrated image The offset:

[0131] .

[0132] Correct the image and corrected image Separately around and Rotation yaw angle Degrees to the route direction, correct the image The pixels in the corrected image can be calculated by the translation transformation matrix The image coordinates in .

[0133] Correct the image and corrected image Separately around and Rotation yaw angle degrees to the direction of the route, it is still a standard upright photography, the visual axes are parallel to each other, and the distance from the center of the camera coordinate system to the ground plane is equal, so the corrected image and corrected image The stitching can be considered as a rigid body translation change, correcting the image Any point in the corrected image can be calculated by the following formula The image coordinates in .

[0134] .

[0135] In some embodiments, in step seven, the method of calculating the overlapping area between the registered images includes:

[0136] The above steps will correct the image and corrected image Separately around and Rotation yaw angle After the two registered images are measured to the flight direction, the overlapping areas of the quadrilaterals formed by the four corners of the two registered images are calculated respectively.

[0137] Specifically, after the rotation is completed, the image is calculated based on computer graphics knowledge. A quadrilateral made of four corners With image , a quadrilateral formed by four corners The overlapping area of ​​polygons.

[0138] Add a certain offset δ outward in the overlapping area to obtain the coordinates of the overlapping area.

[0139] The coordinates of the overlapping area are used to create a mask image, and an image containing only the overlapping area is obtained from the registered image, which are recorded as images and .

[0140] In some embodiments, in step eight, the characteristic band image, the bottom information matrix and the water depth matrix are respectively used as inputs for the accelerated robust feature model SURF2 feature point extraction, and then the integral image of the image is calculated, so that the sum of any features in the image area can be quickly calculated, thereby increasing the speed of subsequent operations. SURF uses a box filter (approximately a Gaussian filter) to convolve the integral image to establish a scale space.

[0141] Construct the DOG scale space, including the following:

[0142] The image pyramid is obtained by Gaussian blurring and downsampling the source image. The downsampling calculation formula is as follows:

[0143] .

[0144] Secondly, each layer of the image pyramid is Gaussian blurred using different parameters, and the multiple blurred images obtained constitute the Gaussian pyramid, and the calculation formula is as follows:

[0145] ;

[0146] .

[0147] Then the images of two adjacent Gaussian spaces are subtracted to obtain the DOG image, and its calculation formula is as follows:

[0148] .

[0149] In the formula, represents downsampling, where ; represents the source image, represents the Gaussian scale space after convolution of the original image, represents the convolution operation, represents the scale factor of the Gaussian convolution kernel; represents the Gaussian kernel function, Represents the size of the convolution kernel; represents the scale factor of the adjacent scale space, and is taken .

[0150] Find points of interest in scale space, including the following:

[0151] In scale space, for each pixel, its 26 neighbors at the current scale are compared (8 at the same scale and 9 at the upper and lower scales).

[0152] Delete unstable feature points, including the following:

[0153] First, in order to obtain the precise position of the interest point, it is necessary to perform pixel difference in discrete space, and then fit the three-dimensional quadratic function by analogy to obtain the precise position of the interest point, and finally eliminate the interest points with low contrast.

[0154] At the same time, interest points with large edge responses are discarded.

[0155] Assign the direction of the feature points, including the following:

[0156] Calculate the modulus of each feature point and direction The calculation formula is as follows:

[0157] ;

[0158] .

[0159] In the formula, is the scale space value of the feature point; for each interest point, the main direction is calculated using the Hessian matrix of its surrounding area.

[0160] Generate feature point descriptors, including the following:

[0161] Taking each interest point as the center, we take the surrounding square area and rotate the image area according to the main direction to ensure that the descriptor is invariant to rotation.

[0162] The region is divided into smaller sub-regions, and four basic information are calculated in each sub-region: the sum of gradients, the sum of absolute values ​​of gradients, and the sum of Haar wavelet responses in the x-direction and the y-direction.

[0163] This information is used to construct a 64-dimensional or 128-dimensional descriptor for each interest point.

[0164] In some embodiments, in step nine, the method for screening matching points includes:

[0165] A query index KD tree method is used to preliminarily search for n feature points in the second image to be stitched that match each feature point in the first image to be stitched, forming a one-to-many hypothetical matching set, wherein the first image to be stitched and the second image to be stitched are two images to be stitched.

[0166] Calculate the spatial distance between each feature point in the first image to be stitched and the feature points in the assumed matching set.

[0167] Find the nearest neighbor spatial distance and the next nearest neighbor spatial distance, and calculate the ratio of the nearest neighbor spatial distance to the next nearest neighbor spatial distance ,when When , the two feature points are matched point pairs, and all matching point pairs are found to obtain a matching point pair array.

[0168] Using the query index KD tree method, select and correct images by spatial distance The closest rectified image to each feature point in The n feature matching points in form a one-to-many hypothetical matching set, and the spatial distance is used as the evaluation index of the feature point similarity. The spatial distance includes the Euclidean distance and pixel coordinate distance The weighted sum of is calculated as follows:

[0169] ;

[0170] ;

[0171] in, and are the feature descriptor vectors of two feature points respectively;

[0172] and is the pixel position coordinate.

[0173] Calculate the spatial distance between each feature point and the hypothesized matching set point , and its calculation formula is as follows:

[0174] .

[0175] in, It is the influencing factor of the pixel coordinate distance item of the matching point.

[0176] The matching screening strategy uses the ratio of the nearest neighbor space distance to the next nearest neighbor space distance, and its calculation formula is as follows:

[0177]

[0178] in is the closest spatial distance, Second close space distance; , is a predetermined ratio, then the two corner points are considered as matching point pairs, thus achieving a rough match of all feature points.

[0179] For feature points with local similar attributes, the screening process of matching points is as follows:

[0180] Calculate the distance between matching point pairs, and divide the maximum and minimum distance values ​​into 10 intervals evenly. The frequency of each interval is , then the frequency of the peak interval is , the corresponding interval is , in the interval [i-1,i+1] The matching point pairs in are the correct matching point pairs, and this matching point pair set is the exact matching point pairs to be found.

[0181] Then, the erroneous feature point pairs are eliminated according to the random sampling consistent RANSAC algorithm, so as to calculate the homography change matrix, and the image to be fused is obtained by multiplying the homography transformation matrix by the image.

[0182] In some embodiments, constructing a transformation model according to the matching point pair array in step 10 includes:

[0183] For water areas, when the total number of feature points of the corresponding characteristic band images is not less than N, and the total number of feature points of the bottom information matrix and the total number of feature points of the water depth matrix are both less than N, the transformation model is generated according to the matching point pair array; otherwise, the distance histogram is used to constrain the global similarity attribute to eliminate the erroneous matching point pairs, update the matching point pair array, and generate the transformation model according to the updated matching point pair array.

[0184] After the transformation relationship between images is determined, the corresponding global projection transformation model can be solved based on the transformation relationship between images. Solving the projection transformation model of the image is to solve the transformation matrix between the corresponding pixels of the image. The method for calculating the image homography matrix is ​​the direct linear transformation method (DLT).

[0185] In image stitching, let the feature matching point sets of the two images to be stitched be recorded as X and X′ respectively. Use homography mapping transformation to fit the relationship between the two, which can be expressed as:

[0186] .

[0187] Where (u, v, 1)T is the coordinate of the feature point in X′, (x, y, 1)T is the coordinate of the feature point in X, and H is the homography matrix between images, which represents the transformation relationship between images.

[0188] The derivation process of the DLT algorithm is as follows: Expand the above formula and divide the first two rows by the third row to get:

[0189] .

[0190] .

[0191] Can be organized as:

[0192] .

[0193] in:

[0194] .

[0195] .

[0196] Finally, the final global projection transformation matrix can be obtained by solving the singular value vector of the least squares solution of the equation.

[0197] In step 11, the method for constructing a multi-resolution fusion model based on the Laplacian pyramid is as follows:

[0198] The fusion of each two images first creates two images and Gaussian Pyramid , , and then create the corresponding 4-layer Laplace pyramid image , .

[0199] Make a mask image of the same size as the image , this mask image represents the fusion position, and then find the Gaussian pyramid of the mask image , which represents the fusion weight of each pixel.

[0200] At each scale or resolution, according to the current scale The Laplace pyramid images of the two images , Add them together to finally get the spliced ​​Laplace pyramid image .

[0201] by The lowest resolution image of is used as the starting image, and the highest resolution stitching result is reconstructed, such as Figure 5 shown.

[0202] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by ordinary technicians in this technical field within the essential scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for stitching UAV hyperspectral images based on shallow water bottom feature matching, characterized in that: The following steps are involved: Step 1: Obtain the drone hyperspectral images to be stitched, and simultaneously obtain the geographic coordinate projection information corresponding to each image and the latitude and longitude corresponding to each image; Step 2, flare data correction step, calculates the radiation value of each UAV hyperspectral image, and uses the radiation value to correct the flare data in the UAV hyperspectral image; Step 3: Invert the bottom information and water depth corresponding to each pixel of each UAV hyperspectral image, and generate a matrix with the same size as the original image; Step 4: geometrically correct the UAV hyperspectral image according to the latitude and longitude information to obtain a corrected image, and add geographic coordinate projection information to the corrected image; Step 5: resample all the images to be stitched to a uniform resolution; Step 6: Select a corrected image sequence, perform rough image registration according to the geographic coordinate projection information of each corrected image, and obtain a registered image of the corrected image and a translation transformation matrix between adjacent registered images; Step 7, calculating the overlapping area between the registered images, and extracting the overlapping area image to obtain the image to be stitched; Step 8, using deep learning technology to automatically identify and extract key features of land and water areas, separate the land area and water area in the image to be stitched, and process the land area and water area respectively to obtain feature points of each land feature band image, as well as feature points of each water feature band image, feature points of the bottom information matrix, and feature points of the water depth; Step nine, respectively determining the total number of feature points of each feature band image of the land area and the total number of feature points of each feature band image of the water area; when the total number is not less than N, establishing a hypothetical matching set of feature points between two adjacent images to be stitched, and using spatial distance as a similarity evaluation index of feature points, screening out matching points whose spatial distances meet the set value, and obtaining a matching point pair array; Step ten, constructing a transformation model according to the matching point pair array; Step 11, calculating the new coordinates of the image according to the transformation model, and splicing the two adjacent images to be spliced; Step 12: After all the images to be stitched are stitched together, the hyperspectral data, water depth data and bottom data are integrated to output the final stitched image.

2. The method for stitching UAV hyperspectral images based on shallow water bottom feature matching according to claim 1 is characterized in that: The step eight comprises: For the land area, the image of each characteristic band is extracted respectively to obtain the land characteristic band image, and the accelerated robust feature model SURF1 is constructed for each land characteristic band image to obtain the feature points of each land characteristic band image; for the water body area, the image of each characteristic band is extracted respectively to obtain the water body characteristic band image, and the bottom sediment information matrix and the water depth matrix of the water body area are extracted, and the accelerated robust feature model SURF2 is constructed for each water body characteristic band image to obtain the feature points of each water body characteristic band image, the feature points of the bottom sediment information matrix and the feature points of the water depth.

3. The method for stitching UAV hyperspectral images based on shallow water bottom feature matching according to claim 2 is characterized in that: For water areas, it also includes judging the total number of feature points of the bottom information matrix and the total number of feature points of the water depth matrix. When the total number of feature points of the water body characteristic band image, the total number of feature points of the bottom information matrix and the total number of feature points of the water depth matrix are all less than N, the transformation model cannot be constructed; if the total number of feature points of the water body characteristic band image is greater than N, and the total number of feature points of the bottom information matrix and the total number of feature points of the water depth matrix are less than N, the transformation model is generated according to the matching point pair array; when none of the above situations are met, the query index KD tree method is used for the feature points of the bottom information matrix and the feature points of the water depth matrix, a hypothetical matching set of feature points between two adjacent images to be spliced ​​is established, and spatial distance is used as a similarity evaluation index of the feature points, matching points whose spatial distance meets the set value are screened out, and a second matching point pair array is obtained, and the distance histogram is used to constrain the global similarity attribute to eliminate mismatched point pairs, the matching point pair array is updated, and a transformation model is generated according to the second matching point pair array.

4. The method for stitching UAV hyperspectral images based on shallow water bottom feature matching according to claim 1 is characterized in that: Step 2 includes: The calculation method of the radiation value L(λ) is: L(λ)=F L (l)*(DN r (L)-DN amb (l)); Where λ is the wavelength, F L (λ) is the calibration coefficient, DN r (λ) is the brightness value of the pixel, DN amb (λ) is noise; Methods for correcting flare data in UAV hyperspectral images using radiation values ​​include: The flare data is identified using a convolutional neural network and classified according to the strength of the flare. The overexposed flare areas are masked, and the remaining flare areas are divided into strong flare areas and weak flare areas. The strong flare areas and weak flare areas are corrected separately using surrounding pixels.

5. The method for stitching UAV hyperspectral images based on shallow water bottom feature matching according to claim 1 is characterized in that: Step 1 also includes obtaining the coordinates of the projection center in the ground photogrammetry coordinate system. The geometric correction in step 4 includes: (41) Establish an image space auxiliary coordinate system parallel to the ground photogrammetry system; (42), the image space auxiliary coordinates (X, Y, Z) of the pixel and the ground photogrammetric coordinates (X A ,Y A ,Z A ) between: Among them, X S ,Y S ,Z S are the coordinates of the pixel projection center in the ground photogrammetry coordinate system, and α is the scale factor.

6. The method for stitching UAV hyperspectral images based on shallow water bottom feature matching according to claim 1 is characterized in that: Adding geographic coordinate projection information to the corrected image includes: The remote sensing image coordinates and geographic coordinates are converted using affine matrix parameters, which have 6 parameters, namely X E ,X pixel ,R γ ,Y E ,Y pixel ,R γ , where X E ,Y E Indicates the geographic projection coordinates of the upper left corner pixel of the corrected image, X pixel ,Y pixel Respectively represent the ground resolution of the corrected image pixels in longitude and latitude, R γ Represents the sine value of the rotation angle of the corrected image; The corrected image is rotated γ degrees clockwise around the center point 0, and then the coordinate projection is calculated according to the corresponding longitude and latitude to obtain the coordinates of point 0 (X OE ,Y OE ), X E ,Y E The calculation formula is: Where w and h represent the width and height of the corrected image respectively; Correct the image by dividing by X E ,Y E The calculation formula for the geographic projection coordinates G(XE′,YE′) of any pixel g(row,col) outside is:

7. The method for stitching UAV hyperspectral images based on shallow water bottom feature matching according to claim 1 is characterized in that: The method for rough image registration in step six is: The geographic coordinates of the upper left corners of the two adjacent corrected images I1 and I2 are G1(X E1 ,Y E1 ) and G2(X E2 ,Y E2 ), taking the corrected image I1 as the reference, the offset of the corrected image I2 relative to the corrected image I1 is calculated: The corrected images I1 and I2 are rotated around O1 and O2 by yaw angles γ to the route direction, respectively. The image coordinates of the pixels in the corrected image I2 in the corrected image I1 can be calculated using the translation transformation matrix.

8. The method for stitching UAV hyperspectral images based on shallow water bottom feature matching according to claim 7 is characterized in that: In step seven, the method for calculating the overlapping area between the registered images includes: Calculate the overlapping area of ​​the quadrilateral formed by the four corners of the two registered images respectively; Add a certain offset δ outward from the overlapping area to obtain the coordinates of the overlapping area; The coordinates of the overlapping area are used to create a mask image, and an image containing only the overlapping area is obtained from the registered image, which are respectively recorded as images I1′ and I1′.

9. The method for stitching UAV hyperspectral images based on shallow water bottom feature matching according to claim 2 is characterized in that: In step nine, the matching point screening method includes: Using a query index KD tree method, preliminarily searching for n feature points in the second image to be stitched that match each feature point in the first image to be stitched, forming a one-to-many hypothetical matching set, wherein the first image to be stitched and the second image to be stitched are two images to be stitched; Calculate the spatial distance between each feature point in the first image to be stitched and the feature points in the assumed matching set; Find out the nearest neighbor spatial distance and the second nearest neighbor spatial distance, calculate the ratio r of the nearest neighbor spatial distance to the second nearest neighbor spatial distance. When r < T, the two feature points are a matching point pair. Find out all the matching point pairs to obtain a matching point pair array.

10. The method for stitching UAV hyperspectral images based on shallow water bottom feature matching according to claim 9 is characterized in that: In step ten, constructing a transformation model according to the matching point pair array includes: For the water body area, when the total number of feature points of each corresponding feature band image is not less than N, and the total number of feature points of the bottom substrate information matrix and the total number of feature points of the water depth matrix are both less than N, generate a transformation model according to the matching point pair array; otherwise, use the distance histogram to constrain the global similarity attribute to eliminate the mis-matched point pairs, update the matching point pair array, and generate a transformation model according to the updated matching point pair array.

Images