A method for bathymetric mapping of nearshore waters based on dual UAV video stitching

Abstract

The present invention discloses a method for surveying and mapping water depth in a nearshore sea area based on dual-UAV video splicing. The method comprises the following steps: completing the calibration of the internal parameters of the cameras of the two UAVs; the cameras carried by the two UAVs face the sea and are distributed along the coastline, shooting videos with a spacing maintained to ensure that the camera fields of view have an intersecting part, and recording the motion characteristics of the waves by shooting videos with optical cameras; then, using the video splicing technology, the videos shot by the two UAVs are spliced; an orthophoto image of each frame is obtained based on the position and attitude information of the UAV cameras; a straight line is fixed in a direction perpendicular to the coastline as a research area, and a time stack diagram of orthophoto corrected images of image key frames of the shot video at the position of the fixed straight line is obtained; and finally, the water depth of the nearshore sea area is obtained by using a cBathy water depth estimation method.

Description

Technical Field

[0001] The present invention relates to the technical field of coastal zone surveying and mapping, and in particular to a method for surveying and mapping the depth of nearshore waters based on dual-UAV video splicing. Background Art

[0002] Traditional methods for bathymetric mapping of nearshore waters rely primarily on physical measuring instruments, such as sonar systems, lidar, synthetic aperture radar, and satellite imagery. However, these instruments require significant manpower and financial resources, making surveying very expensive. Consequently, more cost-effective and efficient surveying methods are needed. With the development of photogrammetry, sensors such as fixed cameras or lidar have emerged to record the motion characteristics of waves. However, this method still requires the deployment of cameras and other measuring tools on-site, requiring appropriate site selection, and is subject to site restrictions. Furthermore, installation and removal are costly.

[0003] Nowadays, as the research on autonomous system drones becomes more mature, the advantages of using this airborne system equipped with a camera are becoming increasingly prominent. It is urgent to propose a low-cost, high-efficiency surveying and mapping method using drones and video processing technology. Summary of the Invention

[0004] The purpose of the present invention is to solve the above-mentioned defects in the prior art and provide a measurement method for estimating the water depth in nearshore waters using images from dual drone-mounted cameras.

[0005] The purpose of the present invention can be achieved by taking the following technical solutions:

[0006] A method for surveying and mapping the depth of coastal waters based on dual UAV video splicing, the method comprising the following steps:

[0007] S1. Complete the internal calibration of the cameras of the two drones. Both drones are equipped with GPS-RTK measurement modules and IMU inertial measurement units to record camera position information and camera attitude information respectively.

[0008] S2. Control two drones to fly along the coastline of the sea area to be measured, keep hovering or move at a constant speed in the same direction with the same spacing and speed to shoot video, ensuring that the fields of view of the two cameras overlap;

[0009] S3, stitching the videos captured by the cameras of the two drones to obtain a panoramic video, and then performing image preprocessing on the panoramic video, wherein the image preprocessing includes image grayscale and image filtering;

[0010] S4. Selecting a sea area to be mapped from the panoramic video after image preprocessing, and generating an orthorectified image of the sea area to be mapped using the recorded camera position information and camera attitude information;

[0011] S5. Using the generated orthorectified image, select a straight line in a direction perpendicular to the coastline, and perform image processing on the straight line to obtain a time stack image of all frames;

[0012] S6. Estimate the water depth information of the corresponding pixel coordinate points in the time stack image using the cBathy water depth estimation method.

[0013] Furthermore, in step S1, the camera intrinsic parameters are calibrated by mathematical tools to obtain the intrinsic parameter matrix. where f x and f y It describes the pixel density of the camera sensor in the horizontal direction of the x-axis and the vertical direction of the y-axis in the image coordinate system. o and y o Represents the pixel offset of the camera's optical axis in the image coordinate system. The camera's intrinsic parameters (abbreviated as intrinsics) describe the inherent properties of the camera itself and determine the shape and size of the 2D image it captures from a 3D scene. Therefore, they are important inputs for image processing and computational geometric transformations.

[0014] Furthermore, the process of step S3 is as follows: define the video shot by the first drone as A, and the video shot by the second drone as B. For each video shot by each drone, divide each frame of the video into m parts according to the width and height of the resolution. Each frame of the video is divided into m 2 grids, and i represents the i-th grid. Let F i (t) represents the homography transformation matrix between the tth frame and the t+1th frame of the ith grid of the video, and the path of the video is defined as the multiplication of the ith grid from the 1st homography transformation matrix to the tth homography transformation matrix, and C i (t) means:

[0015] C i (t) = F i (t)·F i (t-1)···F i (1),1≤t≤T-1,3≤T,1≤i≤m 2

[0016] Where T is the total number of frames of a single video, and C iThe video path represented by (t) reflects the motion trajectory of the camera through the transformation between video frames, and also records the jitter generated during the camera movement, which can effectively describe the motion state of the video. Correspondingly, the video path of each grid i in video A is The video path of each grid i in video B is Set the video path optimization formula Θ(P i )as follows:

[0017] Among them, Ω t represents the range of adjacent frames r of the t-th frame, P i (t) represents the smooth video path of the i-th grid of the video, which is the product of the i-th grid of the video from the first homography transformation matrix to the t-th homography transformation matrix, which is obtained through C i (t) Iterative optimization, P i (r) is also the smooth video path of the i-th grid of the video, which is the multiplication of the i-th grid of the video from the first homography transformation matrix to the r-th homography transformation matrix, which is obtained through C i (r) Iterative optimization is performed, and the smooth video path of each grid i in video A is defined as P i A (t), the smooth video path of each grid i in video B is defined as P i B (t), λ represents the overall weight, which is used to balance ||P i (t)-C i (t)|| and||P i (t)-P i (r)||These two items, w t,r Used to maintain motion discontinuity under rapid displacement or scene change, calculated by Gaussian function G: w t,r =G(||rt||)·G(||C i (r)-C i (t)||)

[0018] Among them, the Gaussian function G is defined as x0 is the input variable of the Gaussian function, μ is the mean of the Gaussian function, σ is the standard deviation of the Gaussian function, and the Gaussian function is used to calculate w t,r Can reflect the relationship between the tth frame and the rth frame, w t,r The larger the value, the closer the motion relationship between the t-th frame and the r-th frame is. t,r The smaller the value, the less close the motion relationship between the tth frame and the rth frame. If the scene changes too quickly, C i (r) and C i The difference in (t) will become larger, w t,r Will become smaller, the video path optimization formula Θ(Pi ) makes the smooth video path focus on maintaining consistency with the original video path, rather than letting the smooth video path focus on fast scene switching, so w is introduced t,r The calculation of C is helpful to solve the situation of fast displacement, rotation and scene switching; the iteration of the video path optimization formula is to calculate C i (t) to update P i (t):

[0019]

[0020] Where ξ is the iteration index, N(i) represents the set of neighboring grids of the i-th grid, j is the j-th grid in the set N(i), and Definition E stable (P) is the video path stabilization formula:

[0021]

[0022] Where P = {P i (t)}, is the smooth video path P corresponding to all grids i i The set of (t), E stable (P) Iteratively optimize the smooth video path P to make the smooth video path close to the original video path, reduce clipping and distortion, and make the smooth video path of each grid smoother. For the smooth video path P of video A and video B A and P B , the video splicing formula is as follows:

[0023] E(P A ,P B ,H)=E stable (P A )+E stable (P B )+β·E stitch (P A ,P B ,H)

[0024]

[0025] Among them, β represents the splicing coefficient, E stitch Indicates the video splicing method, using feature matching method, is the Lth feature point of the tth frame in video A, is the Lth feature point of the tth frame in video B, and the feature point is expressed as a matrix where u fp and v fp The meaning of is the horizontal coordinate and vertical coordinate of the feature point in the pixel coordinate system. and Represent feature points and The grid where it is located, H represents the grid and The passing feature points and The calculated homography matrix;

[0026] Assume that the video path P of the video A shot by the first drone is A Without video path stabilization optimization, the video splicing formula is simplified to:

[0027] E(P A ,P B ,H)=E stable (P B )+β·E stitch (P A ,P B ,H)

[0028] The simplified video splicing formula completes the smoothing of the video path during the iteration process, and enables video B to be spliced ​​with video A. Since only the smoothing of the video path P is required, i B (t) iteration, so that the overall computational complexity is reduced. The smoothed video path P of video B is i B (t) After iterating a predetermined number of times, pass P i B (t)·C i B (t) -1 Each grid in each frame of video B is warped and transformed, and the result is fused frame by frame with video A to obtain a panoramic video. Each frame of the panoramic video is then grayscaled and filtered to remove noise.

[0029] Furthermore, the process of step S4 is as follows:

[0030] Let the coordinate point of the world coordinate system be (x w ,y w ,z w ), the navigation coordinate system used is the north-east coordinate system, the unit is meter, the coordinate of the body coordinate system (x b ,y b ,z b ), the carrier coordinate system used is the front lower right coordinate system, and the coordinate point of the drone camera coordinate system is (x c ,y c ,z c ), the camera coordinate system used is the lower right front coordinate system, the unit is meter, and the process of projecting a point in world space into camera coordinates is:

[0031]

[0032] in The camera's external parameters describe the camera's position, direction, and viewing angle in the 3D scene, and determine the angle from which the camera observes the scene. Therefore, it is an important input for tasks such as 3D reconstruction and pose estimation. 3×3 is the rotation matrix, which is determined by the camera attitude information. Let the roll angle of the camera attitude information be roll, the pitch angle be pitch, and the yaw angle be yaw. Then the calculation formula of the rotation matrix is:

[0033]

[0034] The three-dimensional rotation matrix composed of the posture of the camera coordinate system relative to the body coordinate system, the origin of the camera coordinate system coincides with the origin of the body coordinate system, and the body coordinate system can be rotated around its z-axis by angle θ and then around its y-axis by angle Finally, the rotation angle ψ around its x-axis coincides with the three axes of the camera coordinate system, then

[0035] t 3×1 is the translation matrix, which represents the coordinates of the camera in the world coordinate system. It can be obtained by calculating the camera position and the origin of the world coordinate system. The coordinate point of the camera coordinate system (x c ,y c ,z c ) to pixel plane coordinates (u,v,1) is: Where K is the camera intrinsic parameter matrix obtained in step S1, which is used to confirm the projection properties of the camera.

[0036] The transformation formula of the camera imaging model is:

[0037] Through the above transformation, the points of the world coordinate system are projected onto the two-dimensional pixel plane, and then the image is orthorectified.

[0038] Furthermore, the process of generating an orthorectified image of the selected sea area using the camera position information and the camera attitude information in step S4 is as follows:

[0039] Convert the northeastern coordinates to the local along-coast and perpendicular-coastal coordinate system, setting the local along-coastal direction as the X axis and the cross-coastal direction as the Y axis. Rotate the original northeastern coordinates along the plane of the X and Y axes;

[0040] Define a rectangular area along the coast and perpendicular to the coast, and collect 3D space points in the world coordinate system on the X-axis and Y-axis at equal intervals in the rectangular area, that is, collect 3D space points in the world coordinate system along the coastline and perpendicular to the coastline;

[0041] The pixel coordinates corresponding to each point are obtained through the camera imaging model. After uniformly collecting a series of points in three-dimensional space, the corresponding pixel coordinates are obtained. The pixel coordinates are then sampled, interpolated, and rearranged to form a new image, which is the orthorectified image.

[0042] Furthermore, in step S5, for the orthorectified pixel image, let the total number of pixel images be m0, the number of columns be cols, and the number of rows be rows. Let the value range of col be 1 to cols. Starting with the pixel value in the first column and ending with the pixel value in the cols column, the pixel values ​​in the col column from the first image to the m0th image in the orthophoto are sequentially taken and arranged side by side to form a time stack graph. The width of this time stack graph is the time length, the number of columns is m0, the height is the coastline length, the number of rows is rows, and the total number of time stack graphs is cols. Commonly used bathymetric mapping methods are mostly frequency-domain based. The time stack graph obtained by the above steps is useful for analyzing physical quantities such as wave propagation direction, propagation speed, and propagation frequency. This data structure is also suitable for frequency-domain based methods for water depth estimation.

[0043] Furthermore, the process of step S6 is as follows: According to the dispersion relation, in linear waves, there exists the following relationship among the wave angular frequency ω, the wave number k, the water depth h, and the gravitational acceleration g:

[0044] ω 2 =gktanh(kh)

[0045] The calculation formula for water depth h is derived as follows: In addition, there are identities among the wave velocity c, wave number k, wave angular frequency ω, and wave frequency f: ω=2πf, then the water depth formula can also be written as:

[0046] According to the above formula, only two of the four physical quantities, wave velocity c, wave number k, wave angular frequency ω, and wave frequency f, need to be known to derive an estimated value of the water depth. Then, based on the time stack diagram obtained in step S5, the cBathy water depth estimation method is used to estimate the wave frequency f and wave number k, thereby calculating the water depth h. The cBathy water depth estimation method is a two-dimensional frequency domain estimation method that mainly includes signal processing and fault-tolerant processing. Due to the presence of a large number of nonlinear fitting optimization steps and feedback links, multiple factors are weighted for consideration when estimating parameters, and a fault-tolerant verification operation is used on the final parameters, so the algorithm is more robust.

[0047] The present invention has the following advantages and effects compared to the prior art:

[0048] 1. This invention uses dual drones to record video data while moving around to perform nearshore bathymetry mapping. The scope of the surveyed sea area can be continuously updated by the drones moving along the coastline, while simultaneously recording the already mapped area to expand the surveyed area. Existing methods use a single drone hovering in mid-air or using a fixed-point camera for mapping. In comparison, this invention offers a wider mapping range and higher efficiency.

[0049] 2. When a single drone camera records video, the far-flung coastal waters on the left and right sides are far away from the camera and the angle of view is poor, resulting in information loss that affects the mapping results. By introducing two drone cameras, each camera captures the coastal waters farther from the previous drone, ensuring the overall video information is as complete as possible.

[0050] 3. When stitching videos recorded by drones, the present invention uses a method of bundled video paths and joint path optimization. Instead of using a traditional image stitching algorithm that stitches corresponding frames of each video with the same timestamp, the present invention combines the video paths for joint optimization to eliminate jitter and ghosting generated by image stitching. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] The drawings described herein are used to provide a further understanding of the present invention and constitute a part of this application. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the drawings:

[0052] Figure 1 This is a flow chart of a method for bathymetric mapping of nearshore waters based on dual UAV video splicing disclosed in an embodiment of the present invention;

[0053] Figure 2 This is a frame of effect image captured after video splicing in a method for coastal water depth mapping based on dual UAV video splicing disclosed in an embodiment of the present invention;

[0054] Figure 3 It is a sea area range map (200 meters × 200 meters) selected in a nearshore water depth mapping method based on dual UAV video splicing disclosed in an embodiment of the present invention;

[0055] Figure 4 It is an orthorectified map in a method for nearshore water depth mapping based on dual UAV video stitching disclosed in an embodiment of the present invention;

[0056] Figure 5 A time stack diagram in a method for bathymetric mapping of nearshore waters based on dual UAV video stitching disclosed in an embodiment of the present invention;

[0057] Figure 6 A diagram showing the depth estimation results of a nearshore waters depth mapping method based on dual UAV video stitching in an embodiment of the present invention. DETAILED DESCRIPTION

[0058] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings 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 making creative efforts shall fall within the scope of protection of the present invention.

[0059] Example 1

[0060] like Figure 1 As shown, the present embodiment discloses a method for depth mapping of coastal waters based on dual UAV video stitching, comprising the following steps: completing the internal parameter calibration of the cameras of the two UAVs, operating the two UAVs to shoot videos while maintaining a certain distance, performing panoramic video stitching on the videos shot by the two UAVs, obtaining an orthorectified image through the camera information of each frame, demarcating the research area and obtaining a time stack map of the research area, and estimating the depth of the coastal waters using the cBathy water depth estimation method. Figure 1 Detailed introduction to the implementation process:

[0061] Step S1, calculate the internal parameter matrix of the two drone airborne cameras through mathematical tools such as OpenCV library, the internal parameter matrix where f x and f y It describes the pixel density of the camera sensor in the horizontal direction of the x-axis and the vertical direction of the y-axis in the image coordinate system. o and y oRepresents the pixel offset of the camera optical axis in the image coordinate system. A common method is the Zhang Zhengyou calibration method. For convenience, you can also use Matlab's built-in toolbox to complete the calibration of intrinsic parameters.

[0062] Step S2: Since a single drone is far away from the distant cross-coastal sea area, the waves will appear denser. If there is fog or short wind waves on the sea at this time, the signal-to-noise ratio of the drone camera will be too low. Another purpose of controlling the second drone is to enable it to capture the sea area that the first drone cannot capture well, thereby improving the signal-to-noise ratio of the sea area. In the beach area where the water depth needs to be measured, two drones are operated so that their distribution meets certain requirements, such as the two drones flying along the coastline, the cameras of the two drones facing the sea, the shooting range of the cameras of the two drones must have a public field of view, and the two drones keep hovering or moving at a constant speed in the same direction during the video shooting process. The purpose of leaving a public field of view is to facilitate the subsequent video splicing operation to expand the surveying and mapping perspective. The public field of view of the video shot by the two drones accounts for at least 30% of the field of view of the single drone camera, and the speed of the drone movement does not exceed 1 meter per second.

[0063] Step S3: perform panoramic video stitching on the videos shot by the two drones by bundling video paths and jointly optimizing paths. The specific operations are as follows:

[0064] Define the video shot by the first drone as A and the video shot by the second drone as B. For each video shot by each drone, divide each frame of the video into m parts according to the width and height of the resolution. Each frame of the video is divided into m 2 grids, and i represents the i-th grid. Let F i (t) represents the homography transformation matrix between the tth frame and the t+1th frame of the ith grid of the video, and the path of the video is defined as the multiplication of the ith grid from the 1st homography transformation matrix to the tth homography transformation matrix, and C i (t) means:

[0065] C i (t) = F i (t)·F i (t-1)···F i (1),1≤t≤T-1,3≤T,1≤i≤m 2

[0066] Where T is the total number of frames of a single video. Correspondingly, the video path of each grid i in video A is The video path of each grid i in video B is Set the video path optimization formula Θ(P i )as follows:

[0067]

[0068] Among them, Ω t represents the range of adjacent frames r of the t-th frame, which is ±30 frames. i (t) represents the smooth video path of the i-th grid of the video, which is the product of the i-th grid of the video from the first homography transformation matrix to the t-th homography transformation matrix, which is obtained through C i (t) is iteratively optimized, P i (r) is also the smooth video path of the i-th grid of the video, which is the multiplication of the i-th grid of the video from the first homography transformation matrix to the r-th homography transformation matrix, which is obtained through C i (r) is iteratively optimized. The smooth video path of each grid i in video A is defined as P i A (t), the smooth video path of each grid i in video B is defined as P i B (t), λ represents the overall weight, which is used to balance ||P i (t)-C i (t)|| and||P i (t)-P i (r)|| The initial value of these two items is set to 5. If you want the optimized video to be similar to the original video, you can reduce λ. If you want the optimized video to reduce distortion, you can increase λ. t,r Used to maintain motion discontinuity under rapid displacement or scene change, calculated by Gaussian function G: w t,r =G(||rt||)·G(||C i (r)-C i (t)||)

[0069] Among them, the Gaussian function G is defined as x0 is the variable of Gaussian function, μ is the mean value of Gaussian function, which is usually 0, and σ is the standard deviation of Gaussian function, which can be set to 10 as the initial value. The iteration of video path optimization formula is to calculate C i (t) to update P i (t):

[0070]

[0071] Where ξ is the iteration index, N(i) represents the set of neighboring grids of the i-th grid, j is the j-th grid in the set N(i), and Definition E stable (P) is the video path stabilization formula:

[0072]

[0073] Where P = {P i (t)}, is the smooth video path P corresponding to all grids i i (t). By iterating P i (t), the jitter of the video will be effectively removed, which is a prerequisite for achieving stable video splicing. For the smooth video path P of video A and video B A and P B , the video splicing formula is as follows:

[0074] E(P A ,P B ,H)=E stable (P A )+E stable (P B )+β·E stitch (P A ,P B ,H)

[0075]

[0076] Among them, β represents the splicing coefficient, E stitch Indicates the video splicing method, which uses the feature matching method to splice, which is the Lth feature point of the tth frame in video A. is the Lth feature point of the tth frame in video B, and the feature point is expressed as a matrix where u fp and v fp The meaning of is the horizontal coordinate and vertical coordinate of the feature point in the pixel coordinate system. and Represent feature points and The grid where it is located, H represents the grid and The passing feature points and The calculated homography matrix;

[0077] Among them, the video shot by the first drone is selected without video path optimization, and the video shot by the second drone requires path optimization. The video stitching formula is simplified to:

[0078] E(P A ,P B ,H)=E stable (P B )+β·E stitch (P A ,P B ,H)

[0079] Through the simplified video splicing formula, the smooth video path P of video B is i B (t) When the number of iterations is specified in advance, the i B (t)·C i B (t) -1 Each grid in each frame of video B is distorted and transformed, and the obtained results are fused with the frame-by-frame images of video A to obtain a panoramic video, which is then grayscaled and subjected to Gaussian filtering and denoising operations.

[0080] Note that the path optimization for the video here is actually, to some extent, a homography transformation of the image. This transformation will render the camera intrinsic parameters calibrated in step S1 unusable, making it impossible to map 3D world coordinates to 2D pixel coordinates, or even causing large errors in the mapping process. Therefore, it is necessary to exclude the path optimization from one of the drone videos to ensure that it is not affected by the homography transformation and that the camera intrinsic parameters of that drone can be used in subsequent calculations.

[0081] After the video stitching is completed using the joint path optimization method, the entire panoramic video shares an internal reference, which is exactly the internal reference corresponding to the video A shot by the first drone that has not undergone path optimization. The internal reference required for mapping the research range to the two-dimensional plane in the three-dimensional world coordinate system is then satisfied. Then, the orthorectification of the panoramic video will be performed based on the camera position information and camera attitude information of the first drone. If the video B shot by the second drone is selected during the video stitching process without video path optimization, and the video shot by the first drone is selected for video path optimization, then the corresponding video stitching formula becomes:

[0082] E(P A ,P B ,H)=E stable (P A )+β·E stitch (P A ,P B ,H)

[0083] Subsequent panoramic video processing will be based on the data from the second drone.

[0084] Step S4: The transformation matrix for converting the world coordinates to the image coordinate system is as follows:

[0085]

[0086] Where K is the camera intrinsic parameter matrix obtained in step S1, is the external parameter of the camera, R3×3 is the rotation matrix, which is determined by the camera attitude information. Let the roll angle of the camera attitude information be roll, the pitch angle be pitch, and the yaw angle be yaw. Then the calculation formula of the rotation matrix is:

[0087]

[0088] The three-dimensional rotation matrix composed of the posture of the camera coordinate system relative to the body coordinate system, the origin of the camera coordinate system coincides with the origin of the body coordinate system, and the body coordinate system can be rotated around its z-axis by angle θ and then around its y-axis by angle Finally, the rotation angle ψ around its x-axis coincides with the three axes of the camera coordinate system, then

[0089]

[0090] Among them, t 3×1 is the translation matrix, representing the camera's coordinates in the world coordinate system. This can be obtained by calculating the camera position and the origin of the world coordinate system. Through the above transformations, points in 3D space can be projected onto the 2D image plane, and then orthorectification can be performed on the image. Based on the camera position and pose information for the current frame, the area requiring orthorectification can be delineated.

[0091] Because the sea area needs to be determined from image space information, orthorectification of the image is necessary. For ease of representation, the NE coordinates are converted to a local along-coastal and across-coastal coordinate system. Let the local along-coastal direction be the X-axis, and the across-coastal direction be the Y-axis. Simply rotate the original NE coordinate system along the X and Y-axis plane by a certain angle. This angle can be calculated using the local compass azimuth, and the specific angle can be adjusted based on local conditions.

[0092] After the coordinate system is transformed, it is only necessary to define a rectangular area along the coastline and cross-coast directions.

[0093] After the rectangular area is selected, the final orthorectified image needs to be generated in a simple way: three-dimensional space points are collected at equal intervals on the X and Y axes in the previously defined rectangular area, that is, along the coastline and across the coast.

[0094] Through the above method, the three-dimensional coordinates of the points in the area of ​​interest are now fully determined. The camera imaging model can then be used to determine the pixel coordinates corresponding to each point. After uniformly sampling a series of points in three-dimensional space, their corresponding pixel coordinates are also obtained. These pixel coordinates are then interpolated, using either bicubic or bilinear interpolation. These coordinates are then rearranged to form a new image, which is the orthorectified image. It is also important to note that the interval between distance samples in three-dimensional space is generally around 0.5 meters, and generally does not exceed 1 meter. Otherwise, the orthorectified image will lose significant information, affecting the accuracy of the water depth estimation. The sampling interval can be determined by the focal length, resolution, and height of the drone camera above the sea surface. Therefore, the orthorectified image, after resampling the pixel coordinates, contains distance information between pixels.

[0095] Step S5: After the panoramic video has been orthorectified, a straight line is selected perpendicular to the coastline, i.e., across the coastline. This line has fixed coordinates in one dimension of three-dimensional space. This allows a unique solution to be obtained when inversely transforming from image space to three-dimensional space, thereby measuring the three-dimensional position information in image space. A corrected image of the panoramic video is then generated using the selected area and processed to produce a time-stacked image of all frames in that area.

[0096] Since all orthorectified images are available, to facilitate conversion to the data structure required for time series correlation analysis, a column of pixel values ​​(represented in three-dimensional space as the cross-coast direction perpendicular to the coastline) is selected from all orthorectified images to form a time stack image. The specific operation is as follows: for the orthorectified pixel images, let the total number of pixel images be m0, the number of columns be cols, and the number of rows be rows. Let the value of col range from 1 to cols. Starting with the pixel values ​​in the first column and ending with the pixel values ​​in the cols column, the pixel values ​​in the col column from the first image to the m0th image in the orthorectified image are sequentially selected and arranged side by side to form the time stack image. The width of this time stack image is the time length, and the height is proportional to the cross-coast length.

[0097] Step S6: Based on the time stack diagram obtained above, the water depth can be estimated using the cBathy water depth estimation method:

[0098] According to the dispersion relation, in linear waves, there is the following relationship between the wave angular frequency ω, wave number k, water depth h and gravitational acceleration g: ω 2 =gktanh(kh)

[0099] The calculation formula for water depth h is derived as follows:

[0100] In addition, there are identities among the wave velocity c, wave number k, wave angular frequency ω, and wave frequency f: ω=2πf

[0101] Then the water depth formula is further written as:

[0102] Based on the time stack diagram obtained in the previous step, the cBathy water depth estimation method is used to calculate the wave frequency f and wave number k to estimate the water depth h. The schematic diagram of the water depth estimation result is as follows: Figure 6 The cBathy water depth estimation method was published in 2013 in the Journal of Geophysical Research: Oceans, Volume 118, Pages 2595-2609, with the title: cBathy: A robust algorithm for estimating nearshore bathymetry.

[0103] Example 2

[0104] This embodiment further discloses a method for bathymetric mapping of nearshore waters based on dual UAV video stitching, including the following steps:

[0105] S1. Prepare two drones. In order to obtain the position and attitude information of the camera respectively, the two drones need to carry high-precision GPS-RTK devices and IMU inertial measurement units. Use Matlab mathematical tools or OpenCV library to calculate the intrinsic parameter matrix of the drone's airborne camera. The intrinsic parameter matrix where f x and f y It describes the pixel density of the camera sensor in the horizontal direction of the x-axis and the vertical direction of the y-axis in the image coordinate system. o and y o Indicates the pixel offset of the camera optical axis in the image coordinate system.

[0106] S2. In a beach area where water depth mapping is required, operate drones at an altitude of 50 to 150 meters, ensuring that their distribution meets certain requirements. For example, two drones must be distributed along the coastline, their cameras must face the sea, their camera coverage must be within public view, and both drones must simultaneously hover or move at a constant speed in the same direction while recording video. The video must capture periodic waves propagating toward the coast or beach.

[0107] S3. Video data captured by multiple drones are stitched together by bundling video paths and jointly optimizing paths:

[0108] Define the video shot by the first drone as A, and the video shot by the second drone as B. For each video shot by each drone, divide each frame of the video into m parts according to the width and height of the resolution. Each frame of the video is divided into m 2 grids, and i represents the i-th grid. During the calculation process, if the video resolution is too large, m can be set to 10, and if the video resolution is small, m can be set to 4. Let F i (t) represents the homography transformation matrix between the tth frame and the t+1th frame of the ith grid of the video, and the path of the video is defined as the multiplication of the ith grid from the 1st homography transformation matrix to the tth homography transformation matrix, and C i (t) means:

[0109] C i (t) = F i (t)·F i (t-1)···F i (1),1≤t≤T-1,3≤T,1≤i≤m 2

[0110] Where T is the total number of frames of a single video. Correspondingly, the video path of each grid i in video A is The video path of each grid i in video B is Set the video path optimization formula Θ(P i ), which is defined as:

[0111]

[0112] Among them, Ω t represents the range of adjacent frames r of the t-th frame, P i (t) represents the smooth video path of the i-th grid of the video, which is the product of the i-th grid of the video from the first homography transformation matrix to the t-th homography transformation matrix, which is obtained through C i (t) Iterative optimization, P i The initial value of (t) can be equal to C i (t), P i (r) is also the smooth video path of the i-th grid of the video, which is the multiplication of the i-th grid of the video from the first homography transformation matrix to the r-th homography transformation matrix, which is obtained through C i (r) is iteratively optimized. Correspondingly, the smooth video path of each grid i in video A is defined as P i A (t), the smooth video path of each grid i in video B is defined as P i B (t), λ represents the overall weight, and its initial value can be set to 5, which is used to balance ||P i (t)-Ci (t)|| and||P i (t)-P i (r)||These two items, and w t,r Used to maintain motion discontinuity under rapid displacement or scene change, calculated by Gaussian function G:

[0113] w t,r =G(||rt||)·G(||C i (r)-C i (t)||)

[0114] Among them, the Gaussian function G is defined as x0 is the variable of the Gaussian function, μ is the mean of the Gaussian function, which can be 0, and σ is the standard deviation of the Gaussian function. The iteration of the video path optimization formula is calculated by C i (t) to update P i (t):

[0115]

[0116] Where ξ is the iteration index, N(i) represents the set of neighboring grids of the i-th grid, j is the j-th grid in the set N(i), and Definition E stable (P) is the video path stabilization formula:

[0117]

[0118] Where P = {P i (t)}, is the smooth video path P corresponding to all grids i i (t), for the smooth video path P of video A and video B A and P B , the video splicing formula is as follows:

[0119] E(P A ,P B ,H)=E stable (P A )+E stable (P B )+β·E stitch (P A ,P B ,H)

[0120]

[0121] Among them, β represents the splicing coefficient, which is used to balance E stable and E stitch If you want the spliced ​​video to be stable and shake-free, you can set β to a smaller value so that E stableThe weight of E is larger. If you want a good video splicing effect with less ghosting, you can set β to be larger. stitch The weight of E is larger. If you want a good overall effect, you can take β = 0.01, E stitch Indicates the video splicing method, which uses feature matching method to splice. is the Lth feature point of the tth frame in video A, is the Lth feature point of the tth frame in video B, and the feature point is expressed as a matrix where u fp and v fp The meaning is the horizontal coordinate and vertical coordinate of the feature point in the pixel coordinate system. and They represent the feature points and The grid where it is located. H represents the grid and The passing feature points and The calculated homography matrix;

[0122] Assume that the video path P of the video A shot by the first drone is A Without video path stabilization optimization, the video splicing formula is simplified to:

[0123] E(P A ,P B ,H)=E stable (P B )+β·E stitch (P A ,P B ,H)

[0124] Through the simplified video splicing formula, the smooth video path P of video B is i B (t) When the number of iterations is specified in advance, the i B (t)·C i B (t) -1 Each grid in video B is warped and transformed, and the resulting image is fused frame by frame with video A to create a panoramic video. The panoramic video is then preprocessed, primarily by grayscaling the image and performing Gaussian filtering to remove noise. These operations can be implemented using Matlab image processing functions, calling the OpenCV open-source vision library, or writing source code.

[0125] The reason why a certain video is not optimized for video path stabilization is that video path optimization is essentially a distortion transformation or homography transformation, which will change the distortion degree and transmission angle of the video image, resulting in the unavailability of the camera intrinsic parameter matrix K required for the subsequent calculation of the orthorectified area. The video path P of the video A shot by the first drone is selected. A The video path stabilization optimization is not performed so that the subsequent orthorectification can use the flight data of the first drone, such as the camera position information and camera attitude information.

[0126] In addition, if we simply use the traditional image stitching method and the feature point matching method to stitch two videos frame by frame, the stitching effect will be extremely jittery and ghosting, and the video stitching results cannot be used to generate orthorectified images. Therefore, the introduced bundled video path and joint path optimization method can well solve this problem.

[0127] S4. The transformation matrix for converting the world coordinates to the image coordinate system through the projection relationship from the three-dimensional space to the two-dimensional plane is:

[0128] in is the external parameter of the camera, R 3×3 is the rotation matrix, which is determined by the camera attitude information. Let the roll angle of the camera attitude information be roll, the pitch angle be pitch, and the yaw angle be yaw. Then the calculation formula of the rotation matrix is:

[0129]

[0130] The three-dimensional rotation matrix composed of the posture of the camera coordinate system relative to the body coordinate system, the origin of the camera coordinate system coincides with the origin of the body coordinate system, and the body coordinate system can be rotated around its z-axis by angle θ and then around its y-axis by angle Finally, the rotation angle ψ around its x-axis coincides with the three axes of the camera coordinate system, then

[0131] t 3×1 is the translation matrix, which represents the coordinates of the camera in the world coordinate system. It can be obtained by calculating the camera position and the origin of the world coordinate system. Through a series of transformation matrix operations as above, the points in the three-dimensional world coordinate system can be projected onto the two-dimensional plane pixel coordinate system. Based on the camera position and posture information of the current frame, the water depth estimation area to be studied is selected, and the above operations are used to generate an orthorectified image of this area:

[0132] Since the sea area information needs to be determined through the information in the image space, it is necessary to complete the orthorectification of the image. Figure 3The panoramic image shows an area of ​​40,000 square meters, with 200-meter lengths along the coastline and 200-meter lengths along the cross-coastal direction. For convenience, the northeast coordinates are converted to the local along-coastal and cross-coastal coordinate systems. Let the local along-coastal direction be the X-axis and the cross-coastal direction be the Y-axis. The original northeast coordinate system is simply rotated by a certain angle along the plane of the X and Y axes. This angle can be calculated using the local compass azimuth, and the specific angle can be selected based on local conditions.

[0133] After the coordinate system is transformed, it is only necessary to define a rectangular area along the coastline and cross-coast directions.

[0134] After the rectangular area is selected, the final orthorectified image needs to be generated. A simple approach is to collect 3D points at equal intervals along the X and Y axes, i.e., along the coastline and across the coast, in the previously defined rectangular area.

[0135] Through the above method, the three-dimensional spatial coordinates of the area of ​​interest can now be fully determined, and the pixel coordinates corresponding to each point can be obtained using the camera imaging model. After a series of points are evenly collected in three-dimensional space, their corresponding pixel coordinates are also obtained. The pixel coordinates are then interpolated and sampled, and rearranged to reconstruct a new image, which is the orthorectified image. It is also important to note that the interval between distance samples in three-dimensional space is generally around 0.5 meters, which is generally determined by the resolution of the drone's onboard camera and the aircraft's altitude. Therefore, in the orthophoto image after resampling the pixel coordinates, distance information is present between pixels.

[0136] S5. After the panoramic video has been orthorectified, a straight line is selected perpendicular to the coastline (i.e., across the coast). This line has fixed coordinates in one dimension of three-dimensional space. This allows a unique solution to be obtained when inversely transforming from image space to three-dimensional space, thus allowing the measurement of three-dimensional position information in image space. A corrected image of the panoramic video is then generated using the selected region and processed to produce a time stack image of all frames in that region.

[0137] Since all orthorectified images have been obtained, to facilitate conversion to the data structure required for time series correlation analysis, for the orthorectified pixel images, let the total number of pixel images be m0, the number of columns be cols, and the number of rows be rows. Let the value of col range from 1 to cols. Starting with the pixel values ​​in the first column and ending with the pixel values ​​in the cols column, the pixel values ​​in the col column from the first image to the m0th image in the orthorectified image are sequentially selected and arranged side by side. In other words, a column of pixel values ​​(in three-dimensional space, representing the cross-coast direction perpendicular to the coastline) is selected from all orthorectified images to form a time stack image. The width of this time stack image is the time length (in seconds), and the height is proportional to the cross-coast length.

[0138] S6. The cBathy bathymetry method is a two-dimensional frequency-domain estimation method that primarily involves signal processing and fault-tolerance processing. Due to the numerous nonlinear fitting optimization steps and feedback loops involved, multiple weighted factors are considered when estimating parameters, and fault-tolerance checks are applied to the final parameters. This robustness is key. Based on the time stack plot obtained in the previous step, the cBathy method can be used to estimate the water depth h. The cBathy method was published in 2013 in the Journal of Geophysical Research: Oceans, Volume 118, Pages 2595–2609. The article is titled: "cBathy: A robust algorithm for estimating nearshore bathymetry."

[0139] The above embodiments are preferred implementation modes of the present invention, but the implementation modes of the present invention are not limited to the above embodiments. Any other changes, modifications, substitutions, combinations, and simplifications that do not deviate from the spirit and principles of the present invention should be considered as equivalent replacement methods and are included in the scope of protection of the present invention.

Claims

1. A method for bathymetric mapping of nearshore waters based on dual UAV video splicing, characterized in that: The surveying and mapping method comprises the following steps: S1. Complete the internal calibration of the cameras of the two drones. Both drones are equipped with GPS-RTK measurement modules and IMU inertial measurement units to record camera position information and camera attitude information respectively. S2. Control two drones to fly along the coastline of the sea area to be measured, keep hovering or move at a constant speed in the same direction with the same spacing and speed to shoot video, ensuring that the fields of view of the two cameras overlap; S3, stitching the videos captured by the cameras of the two drones to obtain a panoramic video, and then performing image preprocessing on the panoramic video, wherein the image preprocessing includes image grayscale and image filtering; S4. Selecting a sea area to be mapped from the panoramic video after image preprocessing, and generating an orthorectified image of the sea area using the recorded camera position information and camera attitude information; The process of generating orthorectified images of the selected sea area using camera position information and camera attitude information is as follows: Convert the northeastern coordinates to the local along-coast and perpendicular-coastal coordinate system, assuming the local along-coastal direction as the X-axis and the cross-coastal direction as the Y-axis, and rotate the original northeastern coordinates along the plane of the X and Y axes; Define a rectangular area along the coast and perpendicular to the coast, and collect 3D space points in the world coordinate system on the X-axis and Y-axis at equal intervals in the rectangular area, that is, collect 3D space points in the world coordinate system along the coastline and perpendicular to the coastline; The pixel coordinates corresponding to each point are obtained through the camera imaging model. After uniformly collecting a series of points in three-dimensional space, the corresponding pixel coordinates are obtained. The pixel coordinates are then sampled, interpolated, and rearranged to form a new image. This image is the orthorectified image. S5. Using the generated orthorectified image, select a straight line in the direction perpendicular to the coastline and perform image processing to obtain a time stack image of all frames; wherein, for the orthorectified pixel image, let the total number of pixel images be m0, the number of columns be cols, and the number of rows be rows; let the value range of col be 1 to cols, and start with the pixel value of the first column and end with the pixel value of the cols column, sequentially take out the pixel values ​​of the col column from the first image to the m0th image in the orthophoto image and arrange them side by side to form a time stack image, wherein the width of the time stack image is the time length, the number of columns is m0, the height is the coastline length, the number of rows is rows, and the total number of time stack images is cols; S6. Use the cBathy water depth estimation method to estimate the water depth information of the corresponding pixel coordinate points through the time stack image.

2. The method for bathymetric mapping of coastal waters based on dual UAV video splicing according to claim 1 is characterized in that: In step S1, the internal parameter matrix of the drone camera is obtained using an image processing tool, and the internal parameter matrix where f x and f y It describes the pixel density of the camera sensor in the horizontal x-axis and vertical y-axis directions of the image coordinate system. x0 and y0 represent the pixel offset of the camera optical axis in the image coordinate system.

3. The method for bathymetric mapping of coastal waters based on dual UAV video splicing according to claim 1 is characterized in that: The process of step S3 is as follows: Define the video shot by the first drone as A, and the video shot by the second drone as B. For each video shot by each drone, divide each frame of the video into m parts according to the width and height of the resolution. Each frame of the video is divided into m 2 grids, and use i to represent the i-th grid. Let F i (t) represents the homography transformation matrix between the tth frame and the t+1th frame of the ith grid of the video, and the path of the video is defined as the multiplication of the ith grid from the 1st homography transformation matrix to the tth homography transformation matrix, and C i (t) means: C i (t)=F i (t)·F i (t-1)···F i (1),1≤t≤T-1,3≤T,1≤i≤m 2 Where T is the total number of frames of a single video, and the video path of each grid i in video A is set to The video path of each grid i in video B is Video path optimization formula Θ(P i )as follows: Among them, Ω t represents the range of adjacent frames r of the t-th frame, P i (t) represents the smooth video path of the i-th grid of the video, which is the product of the i-th grid of the video from the first homography transformation matrix to the t-th homography transformation matrix, which is obtained through C i (t) Iterative optimization obtains, P i (r) is also the smooth video path of the i-th grid of the video, which is the multiplication of the i-th grid of the video from the first homography transformation matrix to the r-th homography transformation matrix, which is obtained through C i (r) Iterative optimization is performed to obtain the smooth video path of each grid i in video A, which is defined as P i A (t), the smooth video path of each grid i in video B is defined as P i B (t), λ represents the overall weight, which is used to balance ||P i (t)-C i (t)|| and||P i (t)-P i (r)||These two items, w t,r Used to maintain motion discontinuity under rapid displacement or scene change, calculated by Gaussian function G: w t,r =G(||r-t||)·G(||C i (r)-C i (t)||) Among them, the Gaussian function G is defined as x0 is the independent variable of the Gaussian function, μ is the mean of the Gaussian function, and σ is the standard deviation of the Gaussian function; the iteration of the video path optimization formula is to calculate C i (t) to update P i (t): Where ξ is the iteration index, N(i) represents the set of neighboring grids of the i-th grid, j is the j-th grid in the set N(i), and Definition E stable (P) is the video path stabilization formula: Where P = {P i (t)}, is the smooth video path P corresponding to all grids i i (t), for the smooth video path P of video A and video B A and P B , the video splicing formula is as follows: E(P A ,P B ,H)=E stable (P A )+E stable (P B )+β·E stitch (P A ,P B ,H) Among them, β represents the splicing coefficient, E stitch Indicates the video splicing method, using feature matching method, is the Lth feature point of the tth frame in video A, is the Lth feature point of the tth frame in video B, and the feature point is expressed as a matrix where u fp and v fp The meaning of is the horizontal coordinate and vertical coordinate of the feature point in the pixel coordinate system. and Represent feature points and The grid where it is located, H represents the grid and The passing feature points and The calculated homography matrix; Assume that the video path P of the video A shot by the first drone is A Without video path stabilization optimization, the video splicing formula is simplified to: E(P A ,P B ,H)=E stable (P B )+β·E stitch (P A ,P B ,H) Through the simplified video splicing formula, the smooth video path P of video B is i B (t) After iterating a predetermined number of times, Each grid in each frame of video B is distorted and transformed, and the obtained results are fused frame by frame with video A to obtain a panoramic video. Each frame of the panoramic video is then grayscaled and subjected to image filtering and denoising.

4. The method for bathymetric mapping of coastal waters based on dual UAV video splicing according to claim 2 is characterized in that: The process of step S4 is as follows: Let the coordinate point of the world coordinate system be (x w ,y w ,z w ), the navigation coordinate system used is the north-east coordinate system, the coordinates of the body coordinate system (x b ,y b ,z b ), the carrier coordinate system used is the front lower right coordinate system, and the coordinate point of the drone camera coordinate system is (x c ,y c ,z c ), the camera coordinate system used is the lower right front coordinate system, and the process of projecting a point in world space into camera coordinates is: in is the external parameter of the camera, R 3×3 is the rotation matrix, which is determined by the camera attitude information. Assume that the roll angle of the camera attitude information is roll, the pitch angle is pitch, and the yaw angle is yaw. Then the calculation formula of the rotation matrix is: The three-dimensional rotation matrix composed of the posture of the camera coordinate system relative to the body coordinate system, the origin of the camera coordinate system coincides with the origin of the body coordinate system, and the body coordinate system can be rotated around its z-axis by angle θ and then around its y-axis by angle Finally, the rotation angle ψ around its x-axis coincides with the three axes of the camera coordinate system, then t 3×1 is the translation matrix, which represents the coordinates of the camera in the world coordinate system. It can be obtained by calculating the camera position and the origin of the world coordinate system. The coordinate point of the camera coordinate system (x c ,y c ,z c ) to pixel plane coordinates (u,v,1) is: Where K is the camera intrinsic parameter matrix, which is used to confirm the projection properties of the camera; the transformation formula of the camera imaging model is: Through the above transformation, the points of the world coordinate system are projected onto the two-dimensional pixel plane, and then the image is orthorectified.

5. The method for bathymetric mapping of coastal waters based on dual UAV video splicing according to claim 1 is characterized in that: The process of step S6 is as follows: According to the dispersion relation, in linear waves, there is the following relationship between the wave angular frequency ω, wave number k, water depth h and gravitational acceleration g: ω 2 =gktanh(kh) The calculation formula for water depth h is derived as follows: In addition, there are identities among the wave velocity c, wave number k, wave angular frequency ω, and wave frequency f: ω=2πf, the water depth formula is further written as: According to the above formula, only two of the four physical quantities, namely wave velocity c, wave number k, wave angular frequency ω and wave frequency f, need to be known to derive the estimated value of the water depth. Then, based on the time stack diagram obtained in step S5, the cBathy water depth estimation algorithm is used to complete the estimation of wave frequency f and wave number k, thereby calculating the water depth h.

Images