Method and system for three-dimensional reconstruction of ground surface based on multi-view aerial photography of unmanned aerial vehicle

Abstract

The application discloses a ground surface three-dimensional reconstruction method and system based on a multi-view aerial photography of a UAV, relates to the technical field of ground surface three-dimensional reconstruction, and comprises the following steps: acquiring a sequence of aerial photography images containing a highlight area, and synchronously acquiring light center coordinate data, instantaneous speed data and a global solar incident vector at a corresponding moment; extracting a highlight centroid of the highlight area in the sequence of aerial photography images, and tracking a trajectory of the highlight centroid between continuous frames to generate pixel speed data; the method has the beneficial effects that the generation of point cloud elevation collapse and model cavities can be avoided from a spatial calculation source, and the structural integrity and measurability of a three-dimensional ground surface model of a large-scale high-reflectivity area are ensured.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional surface reconstruction technology, and in particular to a method and system for three-dimensional surface reconstruction based on multi-view aerial photography by unmanned aerial vehicles. Background Technology

[0002] In recent years, with the widespread adoption of unmanned aerial vehicles (UAVs), rapidly constructing large-scale 3D models of the Earth's surface using multi-view aerial images has become a core method for acquiring spatial geographic information. Existing conventional 3D reconstruction processes primarily rely on Structure of Motion (SfM) and Multi-View Stereo (MVS) technologies. Their basic principle is to acquire surface images taken by UAVs from multiple perspectives, find corresponding feature points across images from different viewpoints, establish correspondences based on photometric consistency and parallax constraints, and then recover the 3D morphology of the surface through bundle adjustment and depth map densification. This transforms the surface morphology from a 2D image into a measurable 3D model, thereby providing data support for land surveying and spatial planning.

[0003] Existing 3D reconstruction methods suffer from a fundamental flaw: their feature matching and depth calculation algorithms are entirely based on the physical premise of the "Lambertian (diffuse reflection) assumption," which fails to handle specular highlight wander caused by specular reflection on water surfaces or smooth terrain, leading to model voids and elevation collapse. When drones fly over areas with rivers, lakes, or highly reflective surfaces, sunlight produces specular highlights on these non-Lambertian surfaces, and the pixel positions of these highlight spots shift regularly with the drone's spatial perspective. This violates the principle that "physical features should be relatively stationary." The optical phenomenon of "consistent light intensity from multiple perspectives" can cause existing matching algorithms to collapse, misinterpreting it as drift noise and leaving a large number of holes or distorted flying points in areas such as the water surface. Although existing conventional techniques attempt to extract water surface masks using semantic segmentation and rely on the surrounding shoreline for crude elevation interpolation, this is merely image-level post-processing masking. It not only loses the true geometric elevation of the water surface, but also causes interpolation divergence and model topology collapse when facing large areas of water, and can never obtain the true spatial structure of the area from the underlying optical physics. Summary of the Invention

[0004] In view of the above-mentioned prior art, this application is hereby filed. Embodiments of this application provide a method and system for three-dimensional reconstruction of the earth's surface based on multi-view aerial photography by unmanned aerial vehicles (UAVs). This method avoids point cloud elevation collapse and model voids from the source of spatial computation, ensuring the structural integrity and measurability of three-dimensional earth surface models in large-scale, highly reflective areas.

[0005] According to one aspect of this application, a method for three-dimensional reconstruction of the earth's surface based on multi-view aerial photography by unmanned aerial vehicles is provided, including:

[0006] Acquire aerial image sequences containing highlighted areas, and simultaneously acquire the corresponding optical center coordinates, instantaneous velocity data, and global solar incidence vector;

[0007] Extract the specular centroid of the bright region in the aerial image sequence, and track the trajectory of the specular centroid between consecutive frames to generate pixel velocity data. Combine the optical constraint of the normal bisecting the incident light and the line of sight in the specular reflection law with the global solar incident vector and the optical center coordinate data to construct the time derivative relationship of the line of sight. Substitute the instantaneous velocity data and the pixel velocity data into the time derivative relationship to form a specular motion mapping, and inversely calculate the three-dimensional normal vector data of the bright region.

[0008] Calculate the time dynamic variance of multiple three-dimensional normal vector data in a continuous sequence, remove distorted normal vector data that exceeds the preset water body fluctuation threshold, and retain target normal vector data that is less than or equal to the water body fluctuation threshold;

[0009] The target normal vector data is combined with the spatial intersecting ray formed by the specular centroid and the optical center coordinate data to fit and generate the spatial plane equation of the bright region.

[0010] The spatial plane equation is injected into the 3D reconstruction engine as a global rigid constraint. When generating 3D point clouds based on the aerial image sequence, the conventional photometric matching iteration is triggered by the image mask corresponding to the bright area. The discrete point cloud data falling into the image mask is orthogonally projected onto the spatial plane equation along the normal direction to generate the reconstructed 3D model of the ground surface.

[0011] According to another aspect of this application, a three-dimensional surface reconstruction system based on multi-view aerial photography by unmanned aerial vehicles is provided, comprising:

[0012] The multi-source data acquisition module is used to acquire aerial image sequences containing highlighted areas, and simultaneously acquire the corresponding optical center coordinates, instantaneous velocity data, and the calculated global solar incidence vector.

[0013] The reverse calculation module is used to extract the specular centroid of the bright area in the aerial image sequence, and track the trajectory of the specular centroid between consecutive frames to generate pixel velocity data. It combines the optical constraint of the normal bisecting the incident light and the line of sight in the specular reflection law with the global solar incident vector and the optical center coordinate data to construct the time derivative relationship of the line of sight. It substitutes the instantaneous velocity data and the pixel velocity data into the time derivative relationship to form a specular motion mapping, and reversely calculates the three-dimensional normal vector data of the bright area.

[0014] The temporal variance verification module is used to calculate the temporal dynamic variance of multiple three-dimensional normal vector data in a continuous sequence, remove distorted normal vector data that exceeds the preset water body fluctuation threshold, and retain target normal vector data that is less than or equal to the water body fluctuation threshold.

[0015] The spatial reference fitting module is used to combine the target normal vector data with the spatial intersecting ray formed by the specular centroid and the optical center coordinate data to fit and generate the spatial plane equation of the bright region.

[0016] The constraint reconstruction module is used to inject the spatial plane equation as a global rigid constraint into the 3D reconstruction engine. When generating 3D point clouds based on the aerial image sequence, it triggers the image mask corresponding to the bright area to block the conventional photometric matching iteration, and orthogonally projects the discrete point cloud data falling into the image mask onto the spatial plane equation along the normal direction to generate the reconstructed 3D surface model.

[0017] According to another aspect of this application, an electronic device is provided, including a memory and a processor, the memory being used to store computer-executable instructions, and the processor being used to execute the computer-executable instructions, which, when executed by the processor, implement the steps of the method described above.

[0018] According to another aspect of this application, a computer-readable storage medium is provided that stores computer-executable instructions thereon, which, when executed by a processor, implement the steps of the method described above.

[0019] Compared with existing technologies, the method and system for three-dimensional reconstruction of the land surface based on multi-view aerial photography of UAVs according to the embodiments of this application, by inversely calculating the three-dimensional normal vector data and deriving the geometric state using the laws of optical reflection, replaces the conventional photometric consistency matching. It can obtain effective physical spatial orientation information in the highlight areas lacking surface texture features, avoiding the failure of matching corresponding points caused by the wandering of highlight spots. By calculating the time dynamic variance of the three-dimensional normal vector data and comparing and verifying it using a preset water wave threshold, it can identify and eliminate interference features caused by irregular rigid body displacement (such as ships or floating objects moving on the water surface), eliminate data distortion caused by unnatural fluctuations, and ensure that the final retained target normal vector data can objectively reflect the true macroscopic geometric orientation of the target area, thereby improving the reliability of the reconstruction base data.

[0020] In the process of generating 3D point clouds, the conventional matching iteration is blocked by triggering an image mask, and the discrete point cloud is orthogonally projected onto the spatial plane equation along the normal direction. This allows the algorithm to perform geometric correction on the distorted area at the bottom coordinate calculation stage. Compared with post-processing methods that rely on interpolation based on the features of the surrounding shoreline, this application can avoid the generation of point cloud elevation collapse and model voids from the source of spatial calculation, ensuring the structural integrity and measurability of the 3D surface model in a large area of ​​high reflectivity. Attached Figure Description

[0021] The above and other objects, features, and advantages of this application will become more apparent from the more detailed description of the embodiments of this application in conjunction with the accompanying drawings. The drawings are provided to further illustrate the embodiments of this application and form part of the specification. They are used together with the embodiments of this application to explain this application and do not constitute a limitation thereof. In the drawings, the same reference numerals generally represent the same components or steps.

[0022] Figure 1 This is a schematic diagram of the overall process of the three-dimensional reconstruction method of the ground surface based on multi-view aerial photography by UAVs according to the present invention.

[0023] Figure 2 This is a schematic diagram of the logical framework of the three-dimensional reconstruction method of the ground surface based on multi-view aerial photography by UAVs according to the present invention.

[0024] Figure 3 This is a schematic diagram illustrating the extended logic of the three-dimensional reconstruction method of the Earth's surface based on multi-view aerial photography by UAVs, as described in this invention. Detailed Implementation

[0025] Hereinafter, exemplary embodiments according to this application will be described in detail with reference to the accompanying drawings. Obviously, the described embodiments are merely some embodiments of this application, and not all embodiments of this application. It should be understood that this application is not limited to the exemplary embodiments described herein.

[0026] Example 1:

[0027] To address the issues of model voids, elevation collapse, and distortion flypoints caused by specular reflection leading to specular wandering in 3D reconstruction from multi-view drone aerial photography, traditional matching algorithms based on the Lambertian diffuse reflection assumption suffer from multi-view drone aerial image sequences containing highlighted areas. Simultaneously, optical center coordinates, instantaneous velocity data, and the calculated global solar incidence vector are acquired. The specular centroid of the highlighted areas is extracted, and its pixel velocity data across consecutive frames is tracked. The temporal derivative relationship of the line of sight is constructed using specular reflection optical constraints. The instantaneous velocity and pixel velocity are substituted into the specular motion mapping to inversely calculate the 3D normal vector data. Subsequently, the temporal dynamic variance of the normal vector is calculated, and distortion data exceeding the water ripple threshold is removed to retain the target normal vector data. This target normal vector data is then combined with spatial intersecting rays to fit an absolute spatial plane equation. Finally, this is injected as a global rigid constraint into the 3D reconstruction engine, blocking conventional photometric matching iterations within the image mask. Discrete point clouds falling into this area are forcibly orthogonally projected onto the plane equation along the normal direction, generating a reconstructed 3D surface model.

[0028] This invention introduces a reverse analytical mechanism of specular motion mapping, directly utilizing the coupling of optical laws and kinematic data in the underlying data processing stage to achieve elevation analysis and geometric constraints on non-Lambertian surfaces. This avoids the risks of topological divergence and elevation distortion caused by interpolation repair based on surrounding shoreline features. This approach not only avoids the generation of three-dimensional structural voids and flying points in water areas and smooth surfaces from the source of spatial computation, effectively ensuring the rigid structural integrity and measurement accuracy of the three-dimensional surface model, but also reduces the dependence on image surface texture features, improving the reconstruction quality and applicability of UAV aerial surveying in complex optical reflection scenarios.

[0029] Reference Figures 1-3 As an embodiment of the present invention, a method for three-dimensional reconstruction of the ground surface based on multi-view aerial photography by unmanned aerial vehicles is provided, including: S1-S5.

[0030] To facilitate understanding, the following example uses a specific aerial photography scenario to illustrate this embodiment. Assume that a multi-rotor drone is used to perform a multi-view aerial photography mission on a lake and its surrounding surface located in the suburbs of a city. The area to be reconstructed includes the lake surface (approximately...). The drone was equipped with a visible light camera with an equivalent focal length of approximately 2400px. It flew at a constant speed in the east-west direction along a preset "S" shaped route, with an average flight altitude of approximately 150m and a horizontal flight speed of approximately 8m / s. The camera's image acquisition frame rate was approximately 30 frames / second, and the resolution of a single frame image was approximately 4000×3000 pixels. The flight time was around noon, when the solar altitude angle was approximately 75°. Under the sunlight, the lake surface showed visible continuous bright spots. The position of the bright spots shifted continuously with the drone's flight position, and the shape of the bright spots was an irregular ellipse.

[0031] Figure 1 The illustration shows a method for three-dimensional reconstruction of the earth's surface based on multi-view aerial photography by a UAV according to an embodiment of this application, specifically including:

[0032] Reference Figure 1 In step S1, an aerial image sequence containing the highlighted area is acquired, and the corresponding optical center coordinate data, instantaneous velocity data, and global solar incidence vector are acquired simultaneously.

[0033] In step S1, a bright region refers to a set of pixels in an aerial image whose brightness is significantly higher than that of the surrounding Lambertian reflection areas due to strong specular reflection of incident sunlight by the ground surface. These include, but are not limited to, still water surfaces, wet asphalt roads, glass curtain walls, polished metal roofs, and puddles. The identification of bright regions can be performed in real-time based on an image brightness threshold during or after shooting. For example, a set of pixels whose single-channel brightness is higher than a preset brightness threshold (e.g., normalized brightness above 0.85) and whose connected region area is greater than a preset minimum area threshold (e.g., more than 100 pixels) can be marked as a bright region. In an alternative implementation, the identification of bright regions can also be based on a semantic segmentation model (e.g., a water / road surface segmentation model based on U-Net, DeepLab series, etc.) combined with a brightness criterion.

[0034] Aerial image sequence refers to a time-series image collection consisting of multi-view images continuously acquired by a UAV at a preset frame rate during continuous flight along a preset route, covering the area to be reconstructed; optical center coordinate data refers to the three-dimensional position of the camera's optical center in the world coordinate system. This data can be acquired in real time by combining a satellite positioning receiver (such as an RTK-GNSS module) and an inertial measurement unit (IMU) on the UAV, and converted to the camera's optical center through pre-calibrated camera extrinsic parameters; instantaneous velocity data includes two types: translational velocity components and attitude angular velocity components. Translational velocity components refer to the linear velocity of the camera's optical center along the three coordinate axes in the world coordinate system, which can be directly output by the GNSS-IMU combination or calculated from the optical center coordinate difference between adjacent moments; attitude angular velocity components refer to the rotational angular velocity of the camera's body coordinate system around the camera's own three coordinate axes, which can be directly output by the IMU gyroscope. The global solar incidence vector is a unit direction vector from the sun to the Earth's surface, expressed in the world coordinate system. It can be calculated based on a solar calendar model (such as the NRELSPA algorithm) combined with the longitude, latitude, altitude of the shooting area and the shooting UTC time. When the shooting time is short (e.g., no more than 30 minutes) and the solar motion is negligible, this vector can be regarded as remaining constant throughout the entire shooting process.

[0035] Time synchronization between data can be achieved through hardware synchronization using the timestamps provided by the IMU and the exposure trigger signals of the camera. This allows each frame of image to retrieve the optical center coordinates and instantaneous velocity data corresponding to its exposure time after acquisition, and to form a spatiotemporally aligned data frame in conjunction with the pre-calculated global solar incidence vector.

[0036] Using the aforementioned aerial photography scenario example, the drone continuously collects approximately 150 frames of aerial images over a certain flight segment; for the k-th frame, the onboard GNSS-IMU outputs the optical center coordinates at that time. m, the translational velocity component of the k-th frame. m / s, attitude angular velocity components are approximately (0.005, -0.008, 0.012) rad / s; global solar incidence vector calculated using the NRELSPA algorithm based on the latitude and longitude of the shooting area (e.g., 115°E, 32°N) and UTC time:

[0037] In this frame, by performing threshold segmentation and connected component analysis on the single-channel brightness, a connected bright region of about 2000 pixels on the lake surface was identified as the bright region of this frame.

[0038] The above scheme prepares spatiotemporally aligned inputs of image data, spatial location data, kinematic data, and optical path data for the subsequent construction of specular motion mapping. After acquiring the above input data, since conventional photometric consistency matching in the bright area fails due to the change in specular spot position with the viewing angle, the method enters the geometric inversion stage for the bright area. The core of this stage is that instead of searching for traditionally named feature points in the bright area, the phenomenon of "the specular itself moving continuously with the change of the UAV's viewing angle" is interpreted as the normal vector information of the corresponding water surface or mirror area. That is, the specular wandering phenomenon, which is originally regarded as interference in the traditional reconstruction framework, is transformed into a usable geometric inversion input.

[0039] return Figure 1 In step S2, the highlight centroid of the bright area in the aerial image sequence is extracted, and the trajectory of the highlight centroid between consecutive frames is tracked to generate pixel velocity data. The optical constraint of the normal bisecting the incident light and the line of sight in the specular reflection law is combined with the global solar incident vector and optical center coordinate data to construct the time derivative relationship of the line of sight.

[0040] In this step, the specular centroid refers to the geometric representative point of the pixel set within the bright region, which can be extracted using methods such as the arithmetic centroid method, the gray-level weighted centroid method, or the centroid method after brightness thresholding. In one implementation, for the bright regions identified in each frame of the image, the pixel coordinates of the specular centroid are calculated using the gray-level weighted centroid method. This involves using the brightness of each pixel within the highlighted area as a weight to perform a weighted average of its pixel coordinates. This makes the centroid position sensitive to the brightness distribution inside the highlight spot but robust to edge jitter. The trajectory tracking of the highlight centroid between consecutive frames can be achieved based on optical flow methods (such as Lucas-Kanade sparse optical flow), template matching methods, or consecutive frame association methods based on brightness peaks. Then, by differentiating the pixel coordinates of the highlight centroid in adjacent frames and dividing by the time interval between adjacent frames, the pixel velocity data of the highlight centroid in the pixel plane can be obtained.

[0041] The specular reflection law, characterized by the optical constraint that "the normal bisects the incident light and the line of sight," means that at any given moment, the line of sight from the camera's optical center to the physical point corresponding to the specular centroid, the direction of solar incidence, and the direction of the surface normal at that physical point are coplanar, and the normal direction is the angle bisector of the incident direction and the line of sight. When the drone's flight causes the optical center position to change, according to this constraint, the line of sight must undergo a corresponding continuous change to maintain specular reflection geometry, thus mathematically establishing a time derivative relationship for the line of sight direction. Substituting the drone's translational and attitude kinematic data into this time derivative relationship constitutes the so-called specular motion mapping, which creates an explicit algebraic relationship between the specular velocity in the pixel plane and the water surface normal, camera motion, solar incidence, and spatial depth. This algebraic relationship can then be reverse-engineered to extract the normal vector.

[0042] Traditional 3D reconstruction frameworks are based on the assumption of diffuse reflection of Lambertian bodies. The processing of bright areas depends on finding corresponding points with consistent luminosity between images from different viewpoints. However, specular reflection areas violate this assumption. Specular spots appear at different physical points under different viewpoints, causing luminosity consistency matching to fail. The specular reflection law itself provides geometric constraints between specular highlights, line of sight, and normals. Combined with the kinematic data that the UAV can acquire, the interference phenomenon of "spectrum wandering" under the traditional framework can be transformed into measurable information about the geometric shape of the bright area.

[0043] This application does not require searching for corresponding feature points in the bright areas, thus avoiding the matching failure problem caused by the failure of photometric consistency in areas such as water surfaces and glass curtain walls that lack texture features; it does not rely on additional hardware (such as polarization cameras, multi-source active lighting, etc.), but only uses image data from the airborne camera and kinematic data output by the GNSS-IMU; the output is three-dimensional normal vector data consistent with the world coordinate system, which can be directly entered into the subsequent spatial plane equation fitting stage;

[0044] In other implementations, the step of geometric inversion of the highlight region can also be replaced by polarization imaging (calculating the normal vector of the highlight region based on a polarization camera), photometric stereo (based on multi-source active illumination), or deep learning-based normal vector prediction methods (directly outputting the normal vector based on a convolutional neural network trained with labeled samples). The aforementioned alternative methods either require additional imaging hardware, a large number of labeled samples as training data, or further transfer training for cross-scene generalization. Compared with the inverse analysis method based on optical physics and kinematics used in this step, they have higher hardware and data dependencies. This step couples the optical reflection law with kinematic data, and can complete the calculation of normal vector data using only the existing image data and kinematic data from UAV aerial photography.

[0045] By substituting instantaneous velocity data and pixel velocity data into the time derivative relationship to form a specular motion map, the three-dimensional normal vector data of the highlight area is calculated in reverse, including:

[0046] The attitude angular velocity component and translation velocity component are extracted from the instantaneous velocity data. Based on the attitude angular velocity component, a spatial rotation transformation matrix between adjacent frames is constructed. The spatial rotation transformation matrix is ​​used to perform reverse rotation compensation on the pixel velocity data to remove the field of view movement error introduced by the UAV attitude deflection and obtain the net pixel velocity data.

[0047] Specifically, constructing the spatial rotation transformation matrix between adjacent frames based on the attitude angular velocity components can be achieved using the Rodrigues rotation formula: [The formula is used to] calculate the time interval between adjacent frames. (s), the attitude angular velocity vector ω (rad / s) and Multiplying them yields the rotation vector between the adjacent frames. (rad), and then convert the rotation vector into a 3×3 orthogonal rotation matrix R according to the Rodrigues formula; in another implementation, the attitude angular velocity can also be integrated into a quaternion, and then the quaternion can be converted into a rotation matrix; using the rotation matrix to perform reverse rotation compensation on the pixel velocity data means: applying the rotation matrix to the direction vector in the normalized camera coordinate system corresponding to the pixel velocity data to reverse cancel the rotation effect, thereby removing the influence of the overall field of view movement caused by the UAV attitude deflection on the specular centroid pixel displacement, and the net pixel velocity data obtained after the above compensation. k is determined solely by the translation of the camera's optical center and the geometric properties of the physical point corresponding to the specular centroid, thus decoupled from the drone's attitude jitter.

[0048] Substituting the net pixel velocity data and translational velocity components into the time derivative relationship, and using the global solar incidence vector and optical center coordinate data as spatial optical path constraint boundaries, a joint specular motion mapping equation is constructed with the three-dimensional normal vector of the highlight region and spatial depth as the two unknown parameters. The specific formula is as follows:

[0049] ;

[0050] in, Let be the joint specular motion mapping residual vector (px / s) of the k-th frame. is the net pixel velocity data (px / s) for the k-th frame, and f is the equivalent focal length of the camera (px). Let be the spatial depth (m) of the highlight region corresponding to the k-th frame. This is the orthogonal projection transformation matrix in the camera coordinate system, with values ​​of diag(1, 1, 0). Let n be the translational velocity component (m / s) of the k-th frame, and n be the three-dimensional normal vector (dimensionless). This is the global solar incident vector (dimensionless). The zero vector (px / s) under the ideal mapping;

[0051] In this formula, the reason for combining the three-dimensional normal vector n with... The reason why the specular velocity on the pixel plane is included as two unknown parameters in the above equation is that it is simultaneously affected by the water surface normal vector, camera translation, the direction of sunlight incidence, and the distance from the physical point corresponding to the centroid of the specular to the camera optical center. This reflects the line-of-sight deflection contribution determined by both the water surface normal and the direction of solar incidence;

[0052] A global algebraic cost function is constructed by extracting the joint specular motion mapping equation for multiple frames within a continuous sliding time window. The specific formula is as follows:

[0053] ;

[0054] in, The global algebraic cost function ( K is the total number of image frames within the continuous sliding time window. for The Euclidean second norm square operator is used to characterize the overall error energy of the residual vector. );

[0055] It should be noted that a continuous sliding time window refers to a time window containing K consecutive images selected from an aerial image sequence with a preset step size. The value of the window length K should balance the stability and timeliness of the solution. In one implementation, K can be between 10 and 30 frames, for example, K=20 frames. Within each sliding time window, the squared L2 norms of the residual vectors of the K joint specular motion mapping equations contained in the sliding time window are summed to obtain the global algebraic cost function corresponding to the window. Based on the cost function minimization solution of multiple frames, noise amplification and solution instability problems that may occur in single-frame solution can be suppressed;

[0056] A regularization penalty constraint representing the inter-frame smoothness of the 3D normal vector is embedded in the global algebraic cost function, and a minimization solution is performed: In each round of iterative computation, the 3D normal vector output in the previous round is fixed, the spatial depth corresponding to each frame is calculated by linear least squares, and the calculated spatial depth is used as a known parameter to substitute back into the joint specular motion mapping equation.

[0057] Specifically, the regularization penalty constraint for the smoothness of the three-dimensional normal vector between frames is used to suppress the drastic jumps in the normal vector between frames due to noise amplification during the solution process. In one implementation, the regularization penalty constraint can be the L2-norm weighted sum of the differences between the normal vectors of adjacent frames within the sliding time window, so that while minimizing the principal term of the cost function, the solution result of the normal vector is constrained to have inter-frame continuity.

[0058] The reason for adopting the alternating iterative solution method of "fixing the normal vector to solve for the spatial depth and fixing the spatial depth to solve for the normal vector" is that the joint specular motion mapping equation is nonlinear with respect to the three-dimensional normal vector n (involving...). The terms are linear; if nonlinear optimization is performed on both parameters simultaneously, the solution dimension is high and it is easy to converge to a local minimum; alternating iteration decomposes the original problem into two low-dimensional subproblems: when the normal vector is fixed, the spatial depth has a closed linear least square solution, and the spatial depth of the highlight region in each frame within the sliding time window can be directly calculated frame by frame; when the spatial depth is fixed, the normal vector is updated by the tangent space optimization method described below; the advantage of this solution strategy is that it reduces the nonlinearity of a single solution and improves the numerical stability. In other implementations, the Levenberg-Marquardt method or the Gauss-Newton method can also be used to perform joint nonlinear optimization on all parameters, but it has a high dependence on the initial value and a large solution dimension;

[0059] Update the 3D normal vector for the current round, including:

[0060] The calculated spatial depth is used as a known parameter to substitute back into the joint specular motion mapping equation. The Euclidean partial derivative of the global algebraic cost function with respect to the 3D normal vector output from the previous round is then calculated to generate an unconstrained gradient vector. The specific formula is as follows:

[0061] ;

[0062] in, For the previous round (the first round) The unconstrained gradient vector corresponding to the three-dimensional normal vector output by the wheel. ), This is the three-dimensional normal vector (dimensionless) output from the previous round. Indicates to Solving for the gradient operator of partial derivatives in three-dimensional space. The global algebraic cost function ( );

[0063] Using the 3D normal vector output from the previous round, a spatial orthogonal projection matrix is ​​constructed. This matrix is ​​then used to project the unconstrained gradient vector onto a geometric tangent plane perpendicular to the 3D normal vector output from the previous round. The radial component in the unconstrained gradient vector that causes deformation of the normal vector's magnitude is removed, generating a tangent space optimization direction vector. The specific formula is as follows:

[0064] ;

[0065] in, For the current round, find the optimal direction vector in the tangent space after stripping the radial component. ), It is a third-order identity matrix. The transpose operator for a vector or matrix;

[0066] The reason for projecting the unconstrained gradient vector onto the geometric tangent plane perpendicular to the previous round's normal vector and performing a rotation update along the unit sphere as described below is that the 3D normal vector, in a physical sense, is a unit vector that only represents direction (with a constant magnitude of 1), and its physical feasible region is the unit sphere centered at the origin, not the entire 3D Euclidean space. If a linear step update is performed directly along the unconstrained gradient vector within Euclidean space, the normal vector will produce a radial component with a magnitude deviating from 1, destroying its unit vector property. Forced normalization of the normal vector is required after each iteration, but this post-processing introduces additional solution drift. This step explicitly constrains the update process within the physical feasible region of the unit sphere: the projection matrix... Its function is to separate the gradient vector from the gradient vector. The radial components in the same direction are retained, but only the tangential components located in the tangential plane are retained, thus obtaining... That is, it corresponds to the reasonable direction of movement on a unit sphere;

[0067] The spherical step radian is calculated by combining the tangent space optimization direction vector with a preset iteration step size parameter. The specific formula is as follows:

[0068] ;

[0069] in, The spherical step radius for the current round. The preset iteration step size parameter ( ), for The modulus length;

[0070] It should be noted that the preset iteration step size parameter Used to optimize the tangent space direction vector Modulus length ( The linear mapping is expressed as a step radian on a unit sphere. In one implementation, The value of follows these principles: an excessively large value may lead to an excessively large step size in a single rotation and iteration divergence; an excessively small value will reduce the convergence speed. In one specific implementation, It can be initially set to make the first round This corresponds to values ​​of approximately 0.01 to 0.05 radians, and an adaptive adjustment strategy based on the decrease in the cost function is adopted during the iteration process—when the global algebraic cost function... Maintain during a continuous decline between adjacent rounds It remains unchanged, but will increase when the global algebraic cost function becomes higher. Halve;

[0071] The 3D normal vector output in the previous round is rotated along the unit sphere by the spherical step radian to the tangent space optimization direction vector, generating an updated 3D normal vector for the current round with a constant magnitude of one. The specific formula is as follows:

[0072] ;

[0073] in, This is the updated 3D normal vector in the current round. The tangent space optimization direction vector after stripping the radial component in the current round;

[0074] The relative descent gradient of the global algebraic cost function between adjacent iterations is calculated using the following formula:

[0075] ;

[0076] in, The relative descent gradient between adjacent iterations (dimensionless);

[0077] The iteration terminates when the relative descent gradient is less than the preset solution convergence threshold, and the three-dimensional normal vector output in the most recent update is confirmed as the three-dimensional normal vector data of the highlighted area.

[0078] It should be noted that in the above iterative solution process, the three-dimensional normal vector output from the previous round must be fixed in each iteration. For the first iteration, since there is no output result from the previous iteration, an initial three-dimensional normal vector needs to be preset as a substitute for the previous round's output. The initial three-dimensional normal vector can be selected in ways including but not limited to: taking a vertical unit vector (0, 0, 1) perpendicular to the horizontal plane as the initial value to fit the prior knowledge that the natural water surface is approximately horizontal in macroscopic terms; or taking the normalized result of the arithmetic mean vector of the normal vectors calculated by solving each frame within the first sliding time window. Correspondingly, since the relative descent gradient between adjacent iterations depends on the values ​​of the global algebraic cost function of the previous and current rounds, the relative descent gradient is calculated starting from the second iteration. The first iteration does not calculate the relative descent gradient and does not use it as a criterion for terminating the iteration.

[0079] It should be noted that the convergence threshold is used to determine whether the iteration has reached a stable convergence state, thereby deciding whether to terminate the iteration. In one implementation, the convergence threshold can be determined based on the requirements of the normal vector accuracy and the timeliness requirements of the scene to be reconstructed: for scenes with high normal vector accuracy requirements (such as water surface elevation measurement), the convergence threshold can be taken as a smaller value (e.g., ...). For scenarios requiring moderate accuracy of the normal vector and fast solution (such as preliminary mapping), a larger convergence threshold can be used (e.g., ...). In one specific implementation, the convergence threshold is determined by... In addition, to prevent the iteration process from failing to meet the convergence threshold for a long time under ill-condition conditions, a maximum number of iterations (such as 100 iterations) can be set as a forced exit condition.

[0080] Using the aforementioned aerial photography example, the highlight centroids of the lake surface in 150 consecutive frames of images from a certain flight segment were extracted frame by frame using the gray-scale weighted centroid method, and then sorted according to the time interval between adjacent frames. The original pixel velocity data is calculated; the attitude angular velocity is used to construct the rotation matrix R between adjacent frames using the Rodrigues rotation formula, and the original pixel velocity data is then compensated by reverse rotation to obtain the net pixel velocity data of the k-th frame. px / s; Select K=20 frames as the continuous sliding time window, and set the preset iteration step size parameter. Pick Solve for the convergence threshold. The normal vector and spatial depth within the window are jointly solved using the alternating iterative method described above. After approximately 35 iterations, the relative descent gradient falls below the convergence threshold, the iteration terminates, and the target 3D normal vector of the region corresponding to the specular centroid is output. The corresponding water surface area is approximately horizontal and slightly tilted to the south. At the same time, the spatial depth sequence of the corresponding window is output, and its arithmetic mean is approximately 148.7m.

[0081] return Figure 1 In step S3, the time dynamic variance of multiple three-dimensional normal vector data under the continuous sequence is calculated, and the specific formula is as follows:

[0082] ;

[0083] in, The time dynamic variance is dimensionless. This represents the total number of three-dimensional normal vector samples extracted from the continuous sequence (i.e., the evaluation window). This refers to the three-dimensional normal vector data calculated at the w-th time step in a continuous sequence. To evaluate the arithmetic mean vector of the 3D normal vector data within the window;

[0084] Distorted normal vector data exceeding the preset water body fluctuation threshold are removed, while target normal vector data less than or equal to the water body fluctuation threshold are retained.

[0085] In this step, the temporal dynamic variance is used to measure the dispersion of the 3D normal vector data around its arithmetic mean vector within the evaluation window. The value of the total number W of 3D normal vector samples extracted from the continuous sequence should take into account both the sensitivity to short-term fluctuations and the ability to suppress random noise. In one implementation, W can be taken as 5 to 20 frames, for example, W=10 frames.

[0086] The reason this application uses time dynamic variance as a screening criterion is that the water surface has geometrical stationarity on a macroscopic time scale, meaning that its macroscopic normal vector should fluctuate slightly around a stable central direction on a time scale of several seconds. In contrast, local normal vector disturbances caused by floating objects (such as leaves, aquatic plants, and floating garbage) or moving vessels will appear as large jumps deviating from the average direction in the solution results. By calculating the time dynamic variance of the normal vector time series data and setting a threshold for screening, the normal vector distortion caused by the above-mentioned irregular rigid body displacement can be distinguished from the real fluctuations of the water surface. Outlier normal vectors caused by irregular rigid body displacement are eliminated, so that the retained target normal vector data can more objectively reflect the real macroscopic geometric orientation of the target area, thereby providing a reliable source of directional data for subsequent spatial plane equation fitting.

[0087] In other implementations, the screening of the normal vector sequence can alternatively employ methods such as RANSAC-based random consistency sampling and elimination, Kalman filter-based state variable smoothing tracking, or Mahalanobis distance-based statistical hypothesis testing. RANSAC is robust to outliers but is sensitive to the number of samples and iterations; Kalman filtering is suitable for online rolling tracking but requires a pre-established state transition model; Mahalanobis distance requires prior estimation of the covariance matrix. In contrast, the time dynamic variance method used in this step has lower computational complexity, the threshold parameter has a clear physical meaning (i.e., directly corresponds to the macroscopic fluctuation amplitude of the water surface), and does not require a pre-established dynamic model.

[0088] It should be noted that the water body fluctuation threshold refers to the maximum variance value allowed for the normal vector to fluctuate around its arithmetic mean vector within the evaluation window. Its value is related to the hydrological conditions, wind force level, water type, and the scene in which the water surface is located. The methods for determining the water body fluctuation threshold include, but are not limited to: First, based on the Beaufort wind scale, the water surface fluctuation level is obtained according to the measured wind speed during the shooting period, and then mapped to a preset variance range; Second, based on pre-calibration of the normal vector calculation results on a still water surface (such as an indoor pool) and a known small-wave water surface, a mapping table of "wind and wave level - variance threshold" is established, and the value is obtained by looking up the table based on the on-site wind speed during flight operations; Third, based on the normal vector sample calculated for a known calm water area at the beginning of this flight, its variance is statistically analyzed as the baseline variance, and amplified by a preset safety factor (such as 1.5 to 3 times) to serve as the water body fluctuation threshold for this flight; In one specific implementation, for inland lakes under wind conditions of level 2 to 3, the water body fluctuation threshold can be approximately... .

[0089] Using the aforementioned aerial photography example, for the 150 consecutive frames of normal vector data output in step S2, the temporal dynamic variance is calculated window by window with W=10 frames as the evaluation window. In one window, because a small boat appears within the coverage area (the boat's bow reflection is incorrectly interpreted as a mirror reflection from the water surface within that window), its temporal dynamic variance is... Solution is approximately Exceeding the preset water body fluctuation threshold (taken as...) The normal vector corresponding to this evaluation window was identified as distorted normal vector data and discarded; the time dynamic variance of the remaining windows was all within approximately [missing information]. Within the range, which is within the water body fluctuation threshold, the corresponding normal vector is retained as the target normal vector data, which can be used as a reliable estimate of the macroscopic geometric orientation of the highlighted area.

[0090] Through the above scheme, the calculated normal vector data was subjected to time-series verification and anomaly removal. After retaining the target normal vector data that reflects the true macroscopic geometric orientation of the target area, the method enters the stage of constructing the spatial geometric reference. The core of this stage is to combine the target normal vector data representing the direction with the spatial ray representing the position, so that the geometric description of the highlighted area is expanded from "only direction" to a complete plane equation that "has both direction and position", thereby providing a measurable global reference surface for subsequent point cloud-level geometric correction.

[0091] return Figure 1 In step S4, the target normal vector data is combined with the spatial intersecting ray formed by the coordinate data of the specular centroid and optical center to fit and generate the spatial plane equation of the highlight region. The specific formula is as follows:

[0092] ;

[0093] ;

[0094] in, Let be the coordinates (m) of the intersection of the spatial intersecting ray and the plane. This represents the optical center coordinates (m). The final spatial depth (m). This is the camera intrinsic parameter matrix (px). The pixel coordinates (px) of the specular centroid. To preserve the target normal vector data, The three-dimensional coordinates of any unknown point on the equation of the spatial plane. (m), The physical distance is zero scalar (m). This indicates that pixel coordinates are converted to normalized camera coordinates (dimensionless).

[0095] It should be noted that spatial intersecting rays refer to those formed by the optical center coordinates at corresponding moments. Starting from the point where the centroid of the highlight pixel coordinates is located... via camera intrinsic parameter matrix The three-dimensional ray extending along the line of sight obtained by inverse mapping is mathematically an extension of the physical point corresponding to the optical center and the specular centroid. Its extension length is determined by the final spatial depth Z* obtained in step S2, thereby anchoring the intersection coordinates of the spatial intersecting ray and the plane. As a spatial point in the highlighted area, the camera intrinsic parameter matrix The spatial depth Z* can be obtained through pre-calibration methods such as Zhang Zhengyou calibration method and checkerboard calibration method, and remains unchanged during flight operations. The final spatial depth Z* is taken from the spatial depth solution of the sliding time window corresponding to the time when the relative descent gradient is less than the convergence threshold in step S2, for example, the arithmetic mean or median of the spatial depth of the specular region of all frames in that window.

[0096] The reason this application uses a combination of "target normal vector + spatial intersecting ray" to fit the spatial plane equation is that, mathematically, the point-normal plane equation only requires two elements—"one direction (normal vector) + one position (a point on the plane)"—to uniquely determine a three-dimensional plane. The target normal vector data n* output by steps S2 and S3 precisely provides the directional element of the plane, while the spatial intersecting ray constructed from the specular centroid pixel coordinates and optical center coordinates provides the positional element on the plane (i.e., the intersection coordinates of the spatial intersecting ray and the plane). Combining the two allows us to obtain the absolute spatial plane equation for the highlighted region without relying on surrounding texture features or semantic segmentation boundaries.

[0097] The advantages of adopting the above scheme are as follows: First, it yields a spatial plane equation with absolute position information, which can be directly used as a global rigid constraint for subsequent point cloud-level geometric correction; Second, the entire fitting process only relies on the intermediate results output by steps S2 and S3 and the pre-calibrated camera intrinsic parameters, without relying on the elevation of the surrounding shoreline or the texture features of the waterside transition area, thus avoiding the interpolation divergence problem caused by the excessive distance of the shoreline in large-area water scenes.

[0098] In other implementations, the fitting of the spatial plane equation can alternatively employ: first, performing least-squares plane fitting directly on the sparse point cloud reconstructed within the highlighted area—however, the point cloud of the water surface under the traditional reconstruction framework itself exhibits sparse and noisy characteristics, limiting the fitting results; second, performing plane extrapolation based on the elevation of the surrounding shoreline—however, when the water area is wide or the shoreline elevation fluctuates significantly, the extrapolation results are sensitive to the geometric accuracy of the shoreline points. This step directly provides the absolute spatial plane equation by combining the normal vector and the ray, avoiding dependence on the original distorted point cloud or shoreline auxiliary data.

[0099] Using the aforementioned aerial photography scenario example, for the sequence of highlight centroids corresponding to 150 consecutive frames of images of a certain flight segment, extract the pixel coordinates of the highlight centroids at a certain frame within the window that was determined to be the target normal vector data in step S3. px, optical center coordinate data m, final spatial depth Z*≈148.7m; camera intrinsic parameter matrix Substitute and After performing a normalized camera coordinate transformation, the intersection coordinates of the intersecting ray and the plane can be obtained by extending along its direction by Z*. m; combined with the target normal vector data retained by S3 The spatial plane equation of the highlighted area on the lake surface can be obtained from the point normal form:

[0100] (m), the plane equation corresponds to the horizontal reference plane of the entire lake surface, and its position and direction can be used as a global rigid constraint for subsequent point cloud correction.

[0101] Through the above scheme, the spatial plane equation of the highlighted area is obtained as the absolute spatial reference. After the absolute spatial reference has been constructed, the method enters the source intervention stage of the 3D point cloud generation process. The core of this stage is to inject the obtained spatial plane equation as a global rigid constraint into the 3D reconstruction engine, so that the conventional photometric matching iteration in the highlighted area is replaced by the rigid constraint, thereby avoiding the reconstruction engine from generating flying points and holes in this area.

[0102] return Figure 1 In step S5, the spatial plane equation is injected into the 3D reconstruction engine as a global rigid constraint, including:

[0103] Extract the normal vector of the spatial plane equation, and calculate the two mutually perpendicular tangent vectors located within the spatial plane equation, including:

[0104] First, it should be noted that a 3D reconstruction engine refers to an image geometry reconstruction framework that can perform Structure for Motion Restoration (SfM) and Multi-View Stereo (MVS) processes on multi-view image sequences and output 3D point clouds, triangular meshes, or textured 3D surfaces. This includes, but is not limited to, open-source or commercial frameworks such as COLMAP, OpenMVS, Pix4D, and Metashape. "Global rigid constraint" refers to the geometric constraint applied to the 3D coordinate update direction of discrete point clouds falling into highlighted areas during the iterative optimization process of the point cloud in the 3D reconstruction engine. This constraint ensures that the update behavior of these point clouds is confined to directions that do not deviate from the spatial plane equations. Image mask. This refers to a binary mask image with the same resolution as the aerial image. Pixels with a value of 1 correspond to the highlighted areas identified in step S1, and pixels with a value of 0 correspond to the non-highlighted areas. The image mask can be generated synchronously after the highlighted areas are identified in step S1 and is transmitted to the 3D reconstruction engine along with the corresponding aerial frame. It serves as the basis for the engine to determine whether to trigger conventional photometric matching iteration or global rigid constraints. The 3D surface model refers to the final output measurable 3D representation data that covers the area to be reconstructed (including highlighted and non-highlighted areas). Its representation form includes, but is not limited to, dense 3D point clouds, triangular mesh surfaces, or textured 3D surface models.

[0105] Extract the global solar incidence vector. Based on the normal vector and the global solar incidence vector, generate the first tangent vector. Based on the normal vector and the first tangent vector, obtain the second tangent vector. The specific formula is as follows:

[0106] ;

[0107] ;

[0108] in, Let be the first tangent vector. The second tangent vector, is the global solar incident vector (dimensionless).

[0109] The reason this application uses the cross product of the normal vector and the global solar incident vector as the first tangent vector, and then the cross product of the normal vector and the first tangent vector as the second tangent vector, to construct the tangent vector basis, instead of using any orthogonal basis (such as the projection of the x / y axes in the world coordinate system onto the plane), is because the solar incident vector S is closely related to the geometric relationship of specular reflection: the first tangent vector Displacement in the direction perpendicular to the "sun-normal" plane subtended by the normal vector and the solar incident vector will not cause a first-order change in the specular pixel coordinates, and the image residual in the corresponding direction is insensitive to displacement; the second tangent vector Located within the "sun-normal" plane and perpendicular to the normal vector, displacement in its direction has a strong first-order effect on the specular pixel coordinates. The image residual in the corresponding direction is sensitive to displacement. By aligning the tangent vector with the specular reflection geometry, the photometric distortion attenuation coefficient described later can be specifically applied to directions susceptible to specular shift. Suppress the direction that is less affected, while suppressing the direction that is less affected. (Direction) Retain the complete update range;

[0110] When iterating the coordinates of discrete point cloud data falling within the image mask, the 3D coordinate update step size of the discrete point cloud data is decomposed into the normal displacement component along the normal vector and the tangential displacement component along the two tangential vectors.

[0111] Specifically, the step size for updating the 3D coordinates of the discrete point cloud data before intervention will be... Target normal vector data along the equation of the space plane and the aforementioned first tangent vector , second tangent vector This set of mutually orthogonal direction bases is orthogonally decomposed, that is, the values ​​are obtained separately. exist , , The formula for calculating the scalar projection components in the three directions, including the normal displacement component along the normal vector direction, is as follows:

[0112] ;

[0113] in, Let be the normal displacement component (m) along the direction of the target normal vector. This represents the step size (m) of the 3D coordinate update of the discrete point cloud data before intervention. , , This constitutes a set of orthonormal bases in three-dimensional space. Before intervention, the step size for updating the three-dimensional coordinates of discrete point cloud data is... Under this set of orthogonal bases, it can be uniquely represented as the vector synthesis of the above three mutually orthogonal components, that is, it satisfies This allows for the application of different constraints to the normal displacement component and the two tangential displacement components.

[0114] The normal displacement component is set to zero, and the tangential displacement component is superimposed onto the current 3D coordinates of the discrete point cloud data to generate updated 3D coordinates, including:

[0115] Extract the first tangential component corresponding to the first tangential vector and the second tangential component corresponding to the second tangential vector from the tangential displacement components. The specific formula is as follows:

[0116] ;

[0117] ;

[0118] in, The first tangential component (m) is represented. This is the second tangential component (m). This represents the step size (m) of updating the three-dimensional coordinates of the discrete point cloud data before any intervention. Let be the first tangent vector. It is the second tangent vector;

[0119] Based on the second tangential component and a preset photometric distortion attenuation coefficient, a scaled second tangential component is generated. The first tangential component and the scaled second tangential component are then superimposed onto the current 3D coordinates of the discrete point cloud data to generate updated 3D coordinates. The specific formula is as follows:

[0120] ;

[0121] in, The updated 3D coordinates (m). The current three-dimensional coordinates (m) of the discrete point cloud data. This is the photometric distortion attenuation coefficient;

[0122] Specifically, the physical meaning of setting the normal displacement component to zero is as follows: In the highlighted area, most of the point cloud displacement in the normal direction calculated by the 3D reconstruction engine comes from the false photometric matching residual driven by specular walk. If the point cloud coordinates are updated by this displacement component, it will cause the point cloud to diverge or collapse along the normal direction. By forcibly setting this component to zero, the discrete point cloud falling into the image mask will no longer be updated in the direction perpendicular to the spatial plane equation, thus being constrained to the vicinity of the spatial plane equation and avoiding point cloud divergence and hole generation.

[0123] The reason this application adopts the method of "zeroing the normal displacement and processing the tangential displacement separately according to direction" instead of completely zeroing the entire 3D coordinate update step or forcibly projecting all coordinates onto the plane equation is that while completely zeroing can avoid flying points, it will also lose the true response of the point cloud to the surface texture and micro-undulations in the tangential direction; completely forced projection will cause the reconstruction result to lose the ability to reflect the geometric gradation of the transition zone between the water surface and the shoreline. This step decomposes the update step into two parts, normal and tangential, applying strong constraints (zeroing) to the normal direction to ensure rigidity, and applying weak constraints (distinguishing) to the tangential direction. / And on (Attenuation) to retain controllable geometric adjustment space, taking into account both rigid constraints and natural geometric transition.

[0124] The advantages of adopting the above scheme are: First, geometric constraints are applied to the highlighted areas from the iterative update stage of the point cloud coordinates, avoiding post-processing repairs after the final model is generated; Second, since the update space is still reserved in the tangential direction, the point cloud of the highlighted area can achieve a natural geometric transition at the boundary with the non-highlighted area, avoiding obvious splicing discontinuities.

[0125] In other implementations, the point cloud processing for the highlighted areas can alternatively employ the following methods: First, completely delete the point cloud of the highlighted areas based on a semantic mask and fill it in using interpolation based on the elevation of the surrounding shoreline—this method is prone to interpolation divergence in large water areas due to excessive shoreline distance; Second, learn to fill the holes in the highlighted areas using a point cloud hole-filling network based on deep learning—this method relies on a large number of labeled samples and does not guarantee the physical flatness of the water surface. This step, through the directional decomposition and differential processing of the coordinate update step size, ensures that the geometric reconstruction results of the highlighted areas directly derive from the constraints of optical physics and kinematic laws, rather than relying on extrapolation from surrounding data or model prediction.

[0126] The photometric distortion attenuation coefficient β has a value range of β∈[0,1] and is used in the second tangent vector. The tangential displacement of the point cloud is controlled and suppressed in the direction of movement. The value of β follows the following principle: when β is 0, The displacement in the direction is set to zero during the update, which is equivalent to applying a rigid constraint in that direction as well; when β is 1, Displacement in the direction is preserved at its original magnitude, equivalent to no additional constraint being applied in that direction; when β takes an intermediate value between 0 and 1, for The directional displacement is scaled proportionally. The value of β is related to the strength of the effect of highlight shift on pixel displacement within the highlighted area: the stronger the highlight shift (e.g., the lower the sun altitude angle, the stronger the water surface specular reflection), the smaller the value of β should be to enhance suppression; the weaker the highlight shift (e.g., the sun altitude angle is close to 90°, the water surface specular reflection is weak), the larger the value of β can be to retain more natural adjustment space. In one specific implementation, β can be 0.2~0.5; in another implementation, β can also be dynamically adjusted according to the statistical highlights intensity of the image within the evaluation window.

[0127] When generating 3D point clouds based on aerial image sequences, the conventional photometric matching iteration is triggered by the image mask corresponding to the highlighted area. The discrete point cloud data falling within the image mask is orthogonally projected onto the spatial plane equation along the normal direction to generate the reconstructed 3D model of the ground surface.

[0128] Specifically, triggering the image mask blocking of the highlighted area to block the conventional photometric matching iteration means that during the point cloud iterative optimization of the image sequence by the 3D reconstruction engine, for each discrete point cloud to be updated, it is first determined whether it falls into the area with an image mask value of 1 based on the projection position of the point cloud on each aerial image; for discrete point clouds that fall into the area with an image mask value of 0, their coordinate updates are still driven by the conventional photometric consistency matching method of the 3D reconstruction engine; for discrete point clouds that fall into the area with an image mask value of 1 (i.e., the highlighted area), their coordinate updates are no longer driven by the conventional photometric matching residual, but by the global rigid constraint provided by the aforementioned spatial plane equation, thereby blocking the conventional photometric matching iteration to avoid the false photometric matching residuals generated by the specular walk polluting the coordinates of these point clouds;

[0129] Projecting discrete point cloud data falling within an image mask orthogonally onto the spatial plane equation along the normal direction means: translating the discrete point cloud along the normal vector direction of the spatial plane equation by a distance equal to the normal vector of the spatial plane equation, so that it falls onto the spatial plane equation. The specific formula is as follows:

[0130] ;

[0131] in, The three-dimensional coordinates (m) of discrete point cloud data after orthogonal projection along the normal direction onto the spatial plane equation. The coordinates (m) of the discrete point cloud data falling within the image mask are given. Let be the coordinates (m) of the intersection of the spatial intersecting ray and the plane. For the target normal vector data, This represents the normal distance (m) from the discrete point cloud data to the spatial plane equation.

[0132] Using the aforementioned aerial photography scenario example, when the 3D reconstruction engine generates point clouds for the aforementioned flight segment image sequence, it synchronously generates corresponding image masks for each frame of aerial photography based on the highlighted areas identified in step S1, and inputs the lake surface spatial plane equation fitted by step S4 as a global rigid constraint into the engine; for the discrete point cloud of the lake surface area falling within the image mask, the update step size is decomposed along the normal vector, the first tangent vector, and the second tangent vector in the above manner, the normal displacement component is set to zero, and the second tangent displacement component is multiplied by the photometric distortion attenuation coefficient β (taken as 0.3) and then superimposed on the current 3D coordinates together with the first tangent displacement component. In the 3D model of the ground surface output after the above processing, the point cloud of the lake surface area fits near the water surface reference plane fitted by step S4, and the "lake surface collapse and void" or "lake surface floating point" phenomenon commonly seen in the traditional reconstruction framework does not appear, while the shoreline area maintains the original geometric transition shape.

[0133] Through the steps S1 to S5 above, a complete method chain is constructed, from optical physics and kinematic data acquisition, inverse calculation of three-dimensional normal vectors, temporal verification, spatial plane fitting to point cloud-level source intervention, providing rigid structural integrity and measurability for three-dimensional surface reconstruction of large-scale high-reflectivity areas.

[0134] Furthermore, since the intensity of water body fluctuations may vary at different physical locations (deep water areas, shallow water areas, inlets, sheltered bays, windward banks, etc.), uniformly applying a rigid constraint with zero normal displacement may flatten the actual low ripples and local micro-undulations on the water surface at certain locations with strong dynamic disturbances. Therefore, this application further proposes, after retaining target normal vector data less than or equal to the water body fluctuation threshold, such as... Figure 3 As shown, it also includes:

[0135] The spatial deflection angle between adjacent time points is calculated based on the target normal vector data in a continuous time series. The specific formula is as follows:

[0136] ;

[0137] in, This represents the spatial deflection angle between time t and the adjacent previous time. and These are the target normal vector data at time t and the adjacent previous time t-1, respectively;

[0138] It should be noted that the above It relies on the target normal vector data of time t and the adjacent previous time t-1. For the first time in the continuous time series arranged in chronological order (denoted as the starting time), since there is no adjacent previous time before it, its spatial deflection angle cannot be calculated according to the above formula. Therefore, it is stipulated that: the starting time does not participate in the calculation of the spatial deflection angle, and the spatial deflection angle sequence starts to be generated from the first time after the starting time; or, the spatial deflection angle of the starting time is assigned to 0 to indicate that there is no measurable directional deflection at the starting time in the absence of an adjacent previous time reference. The spatial deflection angle of subsequent times is still calculated normally according to the above formula. Both of the above processing methods can ensure the logical integrity of the spatial deflection angle sequence.

[0139] The spatial deflection angle is then combined with the instantaneous velocity data at the corresponding moment to construct a normal perturbation response sequence;

[0140] Wherein, the spatial deflection angle between time t and the adjacent previous time. Reflecting the drastic change in the direction of the target normal vector between adjacent time points, the normal perturbation response sequence refers to the spatial deflection angle at each time point. The time-series data sequence is formed by combining the instantaneous velocity data at the same moment; the purpose of including the instantaneous velocity data in this sequence is that the motion of the specular centroid in the pixel plane is affected by both the actual ripples of the water surface and the camera motion. The instantaneous velocity data is used as a reference quantity and together with the spatial deflection angle to form a response sequence, which makes it easier to distinguish the continuous deflection part caused by the stability of the camera motion from the sudden deflection part caused by the instantaneous disturbance of the water surface in the subsequent first-order time-domain difference.

[0141] The normal perturbation response sequence is subjected to first-order time-domain difference calculation. The sequence intervals where the difference result is greater than the preset fluctuation extreme value threshold are extracted as perturbation intensification segments. The phase change frequency data and spatial amplitude data of each perturbation intensification segment are then extracted. The specific formulas are as follows:

[0142] ;

[0143] in, This is the phase change frequency data (rad / s). This represents the spatial deflection angle between time t and the adjacent previous time. This represents the spatial deflection angle between time t-1 and the adjacent previous time. The inherent clock step size (s) for continuous sampling of the system;

[0144] ;

[0145] in, For spatial amplitude data, This represents the set of time indices for the sequence intervals (i.e., perturbation excitation segments) where the difference result is greater than the fluctuation extreme value threshold.

[0146] Specifically, the inherent clock step size of the system's continuous sampling The frame rate is determined by the camera's image capture rate; for example, when the capture frame rate is 30 frames per second... s, the fluctuation extreme value threshold is a threshold parameter used to determine the segment of severe dynamic disturbance in the normal perturbation response sequence. The principle for its value is: if the value is too small, the continuous normal direction change caused by the normal movement of the camera will be misjudged as the segment of intensified perturbation; if the value is too large, the real segment of dynamic disturbance of the water surface will be missed. In one embodiment, the fluctuation extreme value threshold can be determined based on the first difference statistic (such as three times the standard deviation) on the normal vector sample corresponding to the water body fluctuation threshold range determined in the aforementioned S3 step. In another embodiment, for inland lakes under wind conditions of level 2 to 3, the fluctuation extreme value threshold can be taken as about 0.02 to 0.05 rad / s.

[0147] By combining phase change frequency data, spatial amplitude data, and a preset hydrodynamic damping coefficient, the dynamic perturbation coefficient of the specular highlights at different physical locations within the bright region is generated. The specific formula is as follows:

[0148] ;

[0149] in, This represents the hyperbola dynamics perturbation coefficient. This is the hydrodynamic damping coefficient. Indicates the time intervals within the perturbation excitation segment. The arithmetic mean phase change frequency data, For spatial amplitude data, and These represent the reference frequency (rad / s) and reference amplitude (rad) that represent the system limits, respectively.

[0150] It should be noted that the hydrodynamic damping coefficient μ reflects the ability of the observed water body to absorb and dissipate the energy of normal micro-oscillations caused by external disturbances. Its value is related to the viscosity, density, average water depth, and boundary geometry of the water body. In one embodiment, for inland deep lakes and reservoirs, μ can take a smaller value (e.g., 0.3~0.5), which indicates that the water body has a weak damping ability and a wide range of normal micro-perturbation propagation. For shallow ponds and water bodies with obvious aquatic plants or floating objects, μ can take a larger value (e.g., 0.7~0.9), which indicates that the water body has a strong damping ability and normal micro-perturbation decays rapidly. The specific value of μ can be obtained by pre-calibrating the normal micro-perturbation response sequence on a known water body, or by looking up a preset lookup table based on the water body type and water depth parameters.

[0151] Reference frequency With reference amplitude This refers to the maximum value (rad / s) of the normal frequency change and the maximum value (rad) of the normal angle change that can be reliably resolved physically by the image acquisition and processing flow used in this step. These values ​​are used to... Spatial amplitude data The reference value is normalized to obtain the final specular dynamic perturbation coefficient. It is a dimensionless quantity. The upper limit of Nyquist sampling can be determined by the camera's frame rate (e.g., at a frame rate of 30 frames per second). (Corresponding to an angular frequency of approximately 15Hz). The normal deflection amplitude can be determined by the solution accuracy (such as a multiple of the minimum resolvable deflection angle corresponding to the solution accuracy).

[0152] Dynamic confidence grading of the highlighted region is performed by meshing based on the dynamic perturbation coefficient of the specular area. When the discrete point cloud data falling within the image mask is orthogonally projected onto the spatial plane equation along the normal direction, corresponding projection convergence tolerance and spatial smoothing weights are assigned according to different dynamic confidence levels, including:

[0153] It should be noted that the dynamic confidence grading of the highlight area through gridding refers to: dividing the highlight area into two-dimensional grid cells of a preset size along the spatial plane equation (for example, using the two-dimensional plane after the highlight area is projected onto the spatial plane equation as a reference, dividing it into a 1m×1m square grid), and for each grid cell, based on the dynamic perturbation coefficient of the specular highlights falling into that cell. The value is mapped to several preset levels (such as low disturbance, medium disturbance, and high disturbance levels); in one implementation, it can be based on... The value range is preset with three levels of thresholds—for example... ≤0.3 maps to a low disturbance level (corresponding to a deep, calm water zone), 0.3< ≤0.7 is mapped to a medium disturbance level (corresponding to a light wind and wave area). >0.7 is mapped to a high perturbation level (corresponding to the windward wave zone or water inlet); the smaller the size of the grid cell, the higher the spatial resolution of the confidence level, but the corresponding computational cost also increases. The projection convergence tolerance and spatial smoothing weight corresponding to each level can be preset according to the principle of "the stronger the perturbation, the larger the tolerance, and the lower the weight" - the larger the tolerance, the greater the range of deviation of the point cloud from the spatial plane equation in the normal direction; the lower the weight, the less intensity of the point cloud is included in the reference during the surface fitting and smoothing process.

[0154] Extract the dynamic confidence level corresponding to the current physical location, and construct a three-dimensional elastic buffer space on both sides of the normal of the spatial plane equation, using the allocated projection convergence tolerance as the normal boundary threshold.

[0155] Specifically, using the allocated projection convergence tolerance as the normal boundary threshold means: let the normal boundary threshold... Numerically equal to the projection convergence tolerance assigned to the dynamic confidence level corresponding to the current physical location; on both sides of the normal to the spatial plane equation, with the distance from the normal to the spatial plane equation as + and- Two planes parallel to the equation of the spatial plane serve as boundaries. The plate-like spatial region enclosed by these two boundary planes and symmetric about the equation of the spatial plane constitutes a three-dimensional elastic buffer space. That is, the three-dimensional elastic buffer space can be represented as a set of spatial points that satisfy the following equation:

[0156] ;

[0157] in Let be the three-dimensional coordinates (m) of any point in three-dimensional space. Let be the coordinates (m) of the intersection of the spatial intersecting ray and the plane. For the target normal vector data, Let X be the absolute value of the normal distance from point X to the equation of the space plane;

[0158] A three-dimensional elastic buffer space refers to a plate-shaped three-dimensional spatial region with the spatial plane equation as the neutral surface, the allocated projection convergence tolerance as the unilateral normal thickness, and symmetry about the neutral surface. When the absolute value of the initial normal distance data of a discrete point cloud is less than or equal to the normal boundary threshold, the point cloud is considered to be within the three-dimensional elastic buffer space, and the position of the point cloud itself is considered a reasonable representation of the real microscopic undulations of the water surface, without the need to forcibly project it onto the spatial plane equation. When the absolute value of the initial normal distance data of a discrete point cloud is greater than the normal boundary threshold, the point cloud is considered to be outside the three-dimensional elastic buffer space, and is orthogonally projected onto the spatial plane equation along the normal direction as described below.

[0159] The initial normal distance from the discrete point cloud data falling within the image mask to the spatial plane equation is calculated using the following formula:

[0160] ;

[0161] in, The initial normal distance data (m) is given. The coordinates (m) of the discrete point cloud data falling within the image mask are given. To preserve the target normal vector data, Let be the coordinates (m) of the intersection of the spatial intersecting ray and the plane.

[0162] Compare the initial normal distance data with the normal boundary threshold;

[0163] Discrete point cloud data with initial normal distance greater than the normal boundary threshold are orthogonally projected onto the spatial plane equation along the normal direction. The three-dimensional spatial coordinates of discrete point cloud data with initial normal distance less than or equal to the normal boundary threshold are retained. Surface fitting smoothing is then performed on the discrete point cloud data with retained three-dimensional spatial coordinates based on the assigned spatial smoothing weights.

[0164] It should be noted that surface fitting smoothing can employ weighted moving least squares surface fitting based on spatial smoothing weights. That is, for discrete point cloud data that fall within the same grid cell and its neighboring grid cells and retain three-dimensional spatial coordinates, the assigned spatial smoothing weights are used as weighting coefficients for each point cloud in the fitting equation. A local surface equation (such as a quadratic surface equation or a thin plate spline surface equation) is solved, and the elevation values ​​of each point cloud are updated to the surface. In another implementation, surface fitting smoothing can also employ bilateral filtering based on spatial smoothing weights or grid-based Gaussian weighted averaging. The purpose of surface fitting smoothing is to suppress high-frequency jitter of point clouds caused by solution noise while preserving the point cloud's response to the real microscopic undulations of the water surface. This ensures that the final three-dimensional point cloud in the bright area both fits the macroscopic water surface represented by the spatial plane equation and retains local dynamic characteristics.

[0165] The reason why this application adopts the extended processing method of "highlight dynamic perturbation coefficient → dynamic confidence level → differentiated allocation of projection convergence tolerance and spatial smoothing weight" is that the intensity of water surface fluctuations varies at different physical locations: deep water areas are less affected by wind and waves, the water surface is calmer, and the calculated normal vector is closer to the actual water surface direction; shallow water areas or inlets are more affected by the superposition of wind, waves and water flow, the water surface fluctuates more violently, and the calculated normal vector has a large micro-perturbation. If a uniform intensity of normal projection constraint is applied to the entire bright area, the low ripples and undulations that actually exist on the water surface may be flattened in local locations with strong dynamic perturbation.

[0166] The advantages of adopting the above scheme are as follows: First, it applies strong projection constraints to local areas with weak disturbances to ensure the stability of the rigid structure; second, it applies weaker projection constraints to local areas with strong disturbances (preserving the original position of the point cloud in the three-dimensional elastic buffer space) and supplements it with surface fitting smoothing, so that the final three-dimensional point cloud can still reflect the dynamic fluctuations of the local water surface while maintaining macroscopic plane fit; third, the entire disturbance quantification, classification and allocation process only relies on the existing target normal vector data and does not rely on additional hydrological measurement equipment or manual annotation.

[0167] In other embodiments, the differentiated processing of local dynamic characteristics of the water surface can also be alternatively adopted as follows: First, identify areas of severe local fluctuations based on the difference in water reflectance in multispectral images—this method requires additional multispectral imaging equipment; Second, identify areas of severe local fluctuations based on the optical flow field (optical flow density in the bright area)—this method has high requirements for image quality and optical flow stability; This step directly uses the perturbation quantization results of the already calculated normal vector sequence in the time dimension to implement the classification, which is consistent with the main process of this method in terms of data source, and does not require the introduction of additional sensing equipment or image processing process.

[0168] Using the aforementioned aerial photography example, the bright area on the lake surface was divided into a 1m×1m grid. For each grid cell, the spatial deflection angle, phase change frequency data, and spatial amplitude data were calculated based on the time series of the target normal vector within it. Then, the dynamic perturbation coefficient of the specular highlights was calculated. One of the grid units is located near the deep water area in the center of the lake. ≈0.18, mapped to a low disturbance level, corresponding to a projection convergence tolerance of approximately 0.05m and a spatial smoothing weight of 0.9. Discrete point clouds with an absolute normal distance greater than 0.05m falling within this cell are orthogonally projected onto the spatial plane equation along the normal direction. The three-dimensional spatial coordinates of discrete point clouds with an absolute normal distance less than or equal to 0.05m are preserved and incorporated into the surface fitting smoothing with a weight of 0.9. A certain grid cell located at the water inlet on the east side of the lake. The value ≈0.82 is mapped to a high disturbance level, corresponding to a projection convergence tolerance of approximately 0.20m and a spatial smoothing weight of 0.4. Discrete point clouds within this unit with an absolute normal distance greater than 0.20m are orthogonally projected onto the spatial plane equation along the normal direction. The three-dimensional spatial coordinates of discrete point clouds with an absolute normal distance less than or equal to 0.20m are preserved and incorporated into the surface fitting smoothing with a weight of 0.4. After the above differential processing, in the final obtained three-dimensional surface model, the water surface in the deep water area of ​​the lake center approximates a rigid plane, while the inlet area retains the local low undulations caused by water flow disturbance.

[0169] Example 2:

[0170] This is one embodiment of the present invention, which differs from the previous embodiment in that:

[0171] A 3D surface reconstruction system based on multi-view drone aerial photography includes:

[0172] The multi-source data acquisition module is used to acquire aerial image sequences containing highlighted areas, and simultaneously acquire the corresponding optical center coordinates, instantaneous velocity data, and the calculated global solar incidence vector.

[0173] The reverse calculation module is used to extract the specular centroid of the bright areas in the aerial image sequence and track the trajectory of the specular centroid between consecutive frames to generate pixel velocity data. It combines the optical constraint of the normal bisecting the incident light and the line of sight in the specular reflection law with the global solar incident vector and optical center coordinate data to construct the time derivative relationship of the line of sight. It substitutes the instantaneous velocity data and pixel velocity data into the time derivative relationship to form the specular motion mapping and reverse calculates the three-dimensional normal vector data of the bright areas.

[0174] The temporal variance verification module is used to calculate the temporal dynamic variance of multiple three-dimensional normal vector data in a continuous sequence, remove distorted normal vector data that exceeds the preset water body fluctuation threshold, and retain target normal vector data that is less than or equal to the water body fluctuation threshold.

[0175] The spatial reference fitting module is used to combine the target normal vector data with the spatial cross rays composed of the coordinate data of the specular centroid and optical center to fit and generate the spatial plane equation of the bright region.

[0176] The constraint reconstruction module is used to inject the spatial plane equation as a global rigid constraint into the 3D reconstruction engine. When generating 3D point clouds based on aerial image sequences, it triggers the image mask corresponding to the highlighted area to block the conventional photometric matching iteration, and orthogonally projects the discrete point cloud data falling into the image mask onto the spatial plane equation along the normal direction to generate the reconstructed 3D model of the ground surface.

[0177] Example 3:

[0178] In one embodiment of the present invention, which differs from the previous embodiment, the electronic device includes one or more processors and a memory.

[0179] A processor can be a central processing unit (CPU) or other form of processing unit with data processing and / or instruction execution capabilities, and can control other components in an electronic device to perform desired functions.

[0180] The memory may include one or more computer program products, which may include various forms of computer-readable storage media, such as volatile memory and / or non-volatile memory. Volatile memory may include, for example, random access memory (RAM) and / or cache memory. Non-volatile memory may include, for example, read-only memory (ROM), hard disk, flash memory, etc.

[0181] In one example, the electronic device may also include input devices and output devices, which are interconnected via a bus system and / or other forms of connection mechanisms (not shown). In addition, depending on the specific application, the electronic device may include any other suitable components.

[0182] Example 4:

[0183] Embodiments of this application may also be computer-readable storage media storing computer program instructions thereon, which, when executed by a processor, cause the processor to perform the steps described in the "Embodiment 1" section of this specification according to the various embodiments of this application.

[0184] Computer-readable storage media may take the form of any combination of one or more readable media. A readable medium may be a readable signal medium or a readable storage medium. A readable storage medium may, for example, include, but is not limited to, electrical, magnetic, optical, electromagnetic, infrared, or semiconductor systems, apparatuses, or devices, or any combination thereof. More specific examples of readable storage media (a non-exhaustive list) include: electrical connections having one or more wires, portable disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fibers, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof.

[0185] The basic principles of this application have been described above with reference to specific embodiments. However, it should be noted that the advantages, benefits, and effects mentioned in this application are merely examples and not limitations, and should not be considered as essential features of each embodiment of this application. Furthermore, the specific details disclosed above are for illustrative and facilitative purposes only, and are not limitations. These details do not restrict the application from being implemented using the specific details described above.

[0186] The block diagrams of devices, apparatuses, devices, and systems involved in this application are merely illustrative examples and are not intended to require or imply that they must be connected, arranged, or configured in the manner shown in the block diagrams. As those skilled in the art will recognize, these devices, apparatuses, devices, and systems can be connected, arranged, and configured in any manner. Words such as “comprising,” “including,” “having,” etc., are open-ended terms meaning “including but not limited to,” and are used interchangeably with them. The terms “or” and “and” as used herein refer to the terms “and / or,” and are used interchangeably with them unless the context clearly indicates otherwise. The term “such as” as used herein refers to the phrase “such as but not limited to,” and is used interchangeably with it.

[0187] It should also be noted that in the apparatus, equipment, and methods of this application, the components or steps can be disassembled and / or recombined. These disassemblies and / or recombinations should be considered as equivalent solutions of this application.

[0188] The above description of the disclosed aspects is provided to enable any person skilled in the art to make or use this application. Various modifications to these aspects will be readily apparent to those skilled in the art, and the general principles defined herein can be applied to other aspects without departing from the scope of this application. Therefore, this application is not intended to be limited to the aspects shown herein, but rather to be accorded the widest scope consistent with the principles and novel features disclosed herein.

[0189] The above description has been given for purposes of illustration and description. Furthermore, this description is not intended to limit the embodiments of this application to the forms disclosed herein. Although numerous exemplary aspects and embodiments have been discussed above, those skilled in the art will recognize certain variations, modifications, alterations, additions, and sub-combinations thereof.

Claims

1. A method for three-dimensional reconstruction of the earth's surface based on multi-view aerial photography by unmanned aerial vehicles, characterized in that, include: Acquire aerial image sequences containing highlighted areas, and simultaneously acquire the corresponding optical center coordinates, instantaneous velocity data, and global solar incidence vector; Extract the specular centroid of the bright region in the aerial image sequence, and track the trajectory of the specular centroid between consecutive frames to generate pixel velocity data. Combine the optical constraint of the normal bisecting the incident light and the line of sight in the specular reflection law with the global solar incident vector and the optical center coordinate data to construct the time derivative relationship of the line of sight. Substitute the instantaneous velocity data and the pixel velocity data into the time derivative relationship to form a specular motion mapping, and inversely calculate the three-dimensional normal vector data of the bright region. Calculate the time dynamic variance of multiple three-dimensional normal vector data in a continuous sequence, remove distorted normal vector data that exceeds the preset water body fluctuation threshold, and retain target normal vector data that is less than or equal to the water body fluctuation threshold; The target normal vector data is combined with the spatial intersecting ray formed by the specular centroid and the optical center coordinate data to fit and generate the spatial plane equation of the bright region. The spatial plane equation is injected into the 3D reconstruction engine as a global rigid constraint. When generating 3D point clouds based on the aerial image sequence, the conventional photometric matching iteration is triggered by the image mask corresponding to the bright area. The discrete point cloud data falling into the image mask is orthogonally projected onto the spatial plane equation along the normal direction to generate the reconstructed 3D model of the ground surface.

2. The method for three-dimensional reconstruction of the earth's surface based on multi-view aerial photography by unmanned aerial vehicles according to claim 1, characterized in that, After retaining the target normal vector data that is less than or equal to the water body fluctuation threshold, the following is also included: The spatial deflection angle between adjacent moments is calculated based on the target normal vector data in the continuous time series, and the spatial deflection angle is combined with the instantaneous velocity data at the corresponding moment to construct the normal perturbation response sequence; The normal perturbation response sequence is subjected to first-order time-domain difference calculation. The sequence intervals where the difference result is greater than the preset fluctuation extreme value threshold are extracted as perturbation intensification segments. The phase change frequency data and spatial amplitude data of each perturbation intensification segment are extracted. The phase change frequency data, the spatial amplitude data, and the preset hydrodynamic damping coefficient are combined to generate the specular dynamic disturbance coefficients for different physical locations within the high-brightness region. Based on the dynamic perturbation coefficient of the high-brightness region, a dynamic confidence level is performed on the gridded high-brightness region. When the discrete point cloud data falling in the image mask is orthogonally projected onto the spatial plane equation along the normal direction, the corresponding projection convergence tolerance and spatial smoothing weight are assigned according to different dynamic confidence levels.

3. The method for three-dimensional reconstruction of the earth's surface based on multi-view aerial photography by unmanned aerial vehicles according to claim 2, characterized in that, The allocation of the corresponding projection convergence tolerance and spatial smoothing weights includes: Extract the dynamic confidence level corresponding to the current physical location, and construct a three-dimensional elastic buffer space on both sides of the normal of the spatial plane equation, using the allocated projection convergence tolerance as the normal boundary threshold. Calculate the initial normal distance data from the discrete point cloud data falling within the image mask to the spatial plane equation, and compare the initial normal distance data with the normal boundary threshold; Discrete point cloud data with initial normal distance greater than the normal boundary threshold are orthogonally projected onto the spatial plane equation along the normal direction. The three-dimensional spatial coordinates of discrete point cloud data with initial normal distance less than or equal to the normal boundary threshold are retained. Surface fitting smoothing is then performed on the discrete point cloud data with retained three-dimensional spatial coordinates based on the assigned spatial smoothing weights.

4. The method for three-dimensional reconstruction of the earth's surface based on multi-view aerial photography by unmanned aerial vehicles according to claim 1, characterized in that, The inverse solution of the three-dimensional normal vector data includes: The attitude angular velocity component and translation velocity component are extracted from the instantaneous velocity data. A spatial rotation transformation matrix between adjacent frames is constructed based on the attitude angular velocity component. The spatial rotation transformation matrix is ​​used to perform reverse rotation compensation on the pixel velocity data to remove the field of view movement error introduced by the UAV attitude deflection and obtain the net pixel velocity data. Substitute the net pixel velocity data and the translation velocity component into the time derivative relationship, and use the global solar incidence vector and the optical center coordinate data as the spatial optical path constraint boundary to construct a joint specular motion mapping equation with the three-dimensional normal vector and spatial depth of the bright area as the two unknown parameters. The joint specular motion mapping equations of multiple frames within a continuous sliding time window are extracted to construct a global algebraic cost function. A regularization penalty constraint characterizing the inter-frame smoothness of the three-dimensional normal vector is embedded in the global algebraic cost function, and a minimization solution is performed: In each round of iteration, the three-dimensional normal vector output in the previous round is fixed, and the spatial depth corresponding to each frame is calculated by linear least squares. The calculated spatial depth is then used as a known parameter to substitute back into the joint specular motion mapping equation to update the three-dimensional normal vector of the current round. Calculate the relative descent gradient of the global algebraic cost function between adjacent iterations. Terminate the iteration when the relative descent gradient is less than a preset solution convergence threshold, and confirm the three-dimensional normal vector output in the most recent update as the three-dimensional normal vector data of the highlighted region.

5. The method for three-dimensional reconstruction of the earth's surface based on multi-view aerial photography by unmanned aerial vehicles according to claim 4, characterized in that, Updating the three-dimensional normal vector for the current round includes: The calculated spatial depth is substituted back into the joint specular motion mapping equation as a known parameter, and the Euclidean partial derivative of the global algebraic cost function with respect to the three-dimensional normal vector output in the previous round is solved to generate an unconstrained gradient vector. Using the three-dimensional normal vector output from the previous round, construct a spatial orthogonal projection matrix. Using the spatial orthogonal projection matrix, project the unconstrained gradient vector onto a geometric tangent plane perpendicular to the three-dimensional normal vector output from the previous round. Remove the radial component in the unconstrained gradient vector that causes deformation of the normal vector's magnitude, and generate a tangent space optimization direction vector. The tangent space optimization direction vector is combined with the preset iteration step size parameter to calculate the spherical step radian. The three-dimensional normal vector output in the previous round is rotated along the unit sphere according to the spherical step radian to the tangent space optimization direction vector, thereby generating the updated three-dimensional normal vector of the current round with a constant magnitude of one.

6. The method for three-dimensional reconstruction of the earth's surface based on multi-view aerial photography by unmanned aerial vehicles according to claim 1, characterized in that, The step of injecting the spatial plane equation as a global rigid constraint into the 3D reconstruction engine includes: Extract the normal vector of the equation of the spatial plane, and calculate two mutually perpendicular tangent vectors located within the equation of the spatial plane. When iterating the coordinates of the discrete point cloud data falling within the image mask, the three-dimensional coordinate update step size of the discrete point cloud data is decomposed into the normal displacement component along the normal vector and the tangential displacement component along the two tangential vectors. The normal displacement component is set to zero, and the tangential displacement component is superimposed on the current three-dimensional coordinates of the discrete point cloud data to generate updated three-dimensional coordinates.

7. The method for three-dimensional reconstruction of the earth's surface based on multi-view aerial photography by unmanned aerial vehicles according to claim 6, characterized in that: The calculation of the two tangent vectors includes: Extract the global solar incidence vector, generate a first tangent vector based on the normal vector and the global solar incidence vector, and obtain a second tangent vector based on the normal vector and the first tangent vector; The generation of updated 3D coordinates includes: Extract the first tangential component corresponding to the first tangential vector and the second tangential component corresponding to the second tangential vector from the tangential displacement components; Based on the second tangential component and the preset photometric distortion attenuation coefficient, a scaled second tangential component is generated; The first tangential component and the scaled second tangential component are superimposed on the current three-dimensional coordinates of the discrete point cloud data to generate updated three-dimensional coordinates.

8. A three-dimensional surface reconstruction system based on multi-view aerial photography by unmanned aerial vehicles, characterized in that, include: The multi-source data acquisition module is used to acquire aerial image sequences containing highlighted areas, and simultaneously acquire the corresponding optical center coordinates, instantaneous velocity data, and the calculated global solar incidence vector. The reverse calculation module is used to extract the specular centroid of the bright area in the aerial image sequence, and track the trajectory of the specular centroid between consecutive frames to generate pixel velocity data. It combines the optical constraint of the normal bisecting the incident light and the line of sight in the specular reflection law with the global solar incident vector and the optical center coordinate data to construct the time derivative relationship of the line of sight. It substitutes the instantaneous velocity data and the pixel velocity data into the time derivative relationship to form a specular motion mapping, and reversely calculates the three-dimensional normal vector data of the bright area. The temporal variance verification module is used to calculate the temporal dynamic variance of multiple three-dimensional normal vector data in a continuous sequence, remove distorted normal vector data that exceeds the preset water body fluctuation threshold, and retain target normal vector data that is less than or equal to the water body fluctuation threshold. The spatial reference fitting module is used to combine the target normal vector data with the spatial intersecting ray formed by the specular centroid and the optical center coordinate data to fit and generate the spatial plane equation of the bright region. The constraint reconstruction module is used to inject the spatial plane equation as a global rigid constraint into the 3D reconstruction engine. When generating 3D point clouds based on the aerial image sequence, it triggers the image mask corresponding to the bright area to block the conventional photometric matching iteration, and orthogonally projects the discrete point cloud data falling into the image mask onto the spatial plane equation along the normal direction to generate the reconstructed 3D surface model.

9. An electronic device comprising a memory and a processor, characterized in that: The memory is used to store computer-executable instructions, and the processor is used to execute the computer-executable instructions, which, when executed by the processor, implement the steps of the method as described in any one of claims 1 to 7.

10. A computer-readable storage medium having computer-executable instructions stored thereon, characterized in that: When the computer-executable instructions are executed by a processor, they implement the steps of the method as described in any one of claims 1 to 7.

Images