Photogrammetric space reconstruction method based on multi-source image fusion

Abstract

The present application relates to the field of photogrammetry, and discloses a photogrammetry space reconstruction method based on multi-source image fusion, comprising: extracting inner orientation elements and initial outer orientation elements, establishing a meeting constraint strength model of image beams in object space to calculate constraint weights; monitoring the space change rate of the weights to identify the connection area of the heterogeneous images and extract high-weight image beams as geometric anchor points; applying a differential regularization term to the geometric anchor points to construct a stable geometric skeleton, and using a weight damping factor to constrain low-weight image beams to perform conformal mapping correction, wherein the present application characterizes the constraint strength of the image beams, improves the condition number of the adjustment equation matrix, suppresses numerical oscillation caused by differences between heterogeneous sensors, eliminates geometric steps at the fusion interface, and realizes geometric closure and seamless splicing of multi-source images in a unified framework.

Description

Technical Field

[0001] This invention relates to a photogrammetric spatial reconstruction method based on multi-source image fusion, belonging to the field of photogrammetric technology. Background Technology

[0002] Currently, using bundle adjustment to process image sequences to reconstruct the spatial topological relationships of target objects is a mainstream technical approach. This involves establishing the geometric mapping relationship between image points and object points based on the collinearity equation, and solving for the camera exterior orientation elements and object point spatial coordinates by minimizing the sum of squared residuals.

[0003] With the increasing demand for real-scene digital operation and maintenance, integrating UAV high-altitude imagery with ground close-range imagery has become a core operational mode. Due to the significant differences in focal length, pixel physical size, and shooting distance among different source sensors, the geometric constraint strength of the projection ray beams generated by different image sources on spatial points exhibits asymmetrical characteristics. Existing technologies typically use the assumption of homogeneity of observation values. This approach, which ignores the differences in the physical stiffness of image beams, leads to the degradation of the condition number of the global adjustment method equation matrix, resulting in numerical oscillations at the computational level. Although there are abundant hardware acquisition methods for multi-source ground and air imagery, there are still limitations at the level of back-end cross-scale data collaborative adjustment algorithms. For example, Chinese invention patent CN106874436B discloses a multi-source image fusion imaging system for a 3D police geographic information platform, mentioning ground-air joint modeling and overlay of real-scene 3D and panoramic images. The underlying adjustment algorithm relies on POS assistance, DSM guidance, and a conventional feature matching overall process. To address this issue, conventional approaches mainly rely on increasing the number of ground control points or setting empirical weights based on resolution. However, deploying a large number of control points will prolong the engineering cycle, and relying solely on pixel-level probabilistic weighting cannot touch the essential connection between the sensor's physical intrinsic parameters and the geometric stability of the projection. When processing the interface between high-resolution detail images and low-resolution global images, weight mismatch will cause geometric steps, visual tremors, and discontinuities in spatial topological connections in the reconstruction model.

[0004] Therefore, how to establish a cross-scale adjustment weighting system based on physical stiffness constraints and ensure the topological consistency of the model fusion interface has become the technical problem to be solved by this invention. Summary of the Invention

[0005] To address the problems mentioned in the background art, the technical solution of the present invention is as follows: A photogrammetric spatial reconstruction method based on multi-source image fusion, comprising the following steps: Step 1: Acquire UAV-borne images and close-up ground images, and extract the interior orientation elements and initial exterior orientation elements of each image acquisition unit; Step 2: Analyze the focal length in the interior orientation elements and the photographic distance and optical axis tilt angle corresponding to the initial exterior orientation elements. Based on the proportional relationship between focal length and photographic distance and combined with the cosine component of optical axis tilt angle, establish a geometric constraint strength model that characterizes the intersection stability of image bundles in object space, so as to generate the intersection constraint weights of each image bundle in the spatial intersection topology. Step 3: Monitor the spatial change rate of the intersection constraint weights between adjacent image bundles to determine the weight gradient. Based on the weight gradient, identify the connection area between the UAV-borne image and the ground close-range image, and extract the image bundles in the connection area whose intersection constraint weights are higher than a preset threshold as geometric reference anchor points. Step 4: In the iterative optimization process of global bundle adjustment, a second-order differential regularization term based on the spatial position offset is applied to the exterior orientation element corresponding to the geometric reference anchor point. The second-order differential regularization term is used to limit the correction step size of the geometric reference anchor point in the iteration space, thereby constructing a stable geometric skeleton composed of geometric reference anchor points in the object space. Step 5: Construct a numerical damping factor using intersection constraint weights. Use the numerical damping factor to suppress the numerical oscillations generated by the low-weight image bundle during the iteration process. By constructing a local coordinate transformation operator with a stable geometric skeleton as the reference system, force the low-weight image bundle to perform conformal mapping correction under the constraint of the local coordinate transformation operator, thereby locking the geometric offset degree of freedom of the low-weight image bundle relative to the stable geometric skeleton, so as to eliminate the geometric step of the reconstruction model at the image fusion interface.

[0006] Preferably, the logical rule followed in step 2 for calculating the intersection constraint weights w of the image bundle is as follows: obtain the focal length f, the shooting distance D, the optical axis tilt angle θ, and the pixel size δ of the image acquisition unit, and calculate the intersection constraint weights w according to the following formula: w=(f⋅cosθ) / (1000⋅D⋅δ), where w is the intersection constraint weight, f is the focal length in mm, θ is the optical axis tilt angle in ∘, D is the shooting distance in m, and δ is the pixel size in mm; and in the initialization stage of the stochastic model of global bundle adjustment, the magnitude of the normalized intersection constraint weights w is used as the initial input parameter of the observation weight matrix.

[0007] Preferably, step 3 includes the following sub-steps: Step 3-1, using the existence of shared feature points between image bundles as the basis for adjacency determination, establish a topological connection map between adjacent image bundles in the connection area; Step 3-2, calculate the first-order difference of the intersection constraint weights between each node in the topological connection map to obtain the weight gradient; Step 3-3, when the weight gradient exceeds the preset topological mismatch threshold, determine that there is a risk of numerical oscillation of the coefficient matrix of the adjustment method equation induced by the difference in the magnitude of the physical properties of the sensors in the connection area.

[0008] Preferably, applying a second-order differential regularization term based on spatial position offset in step 4 includes: extracting the collinearity constraint equation formed by the geometric reference anchor points in the object space; adding a second-order Laplace penalty term to the column vector of the exterior orientation elements corresponding to the geometric reference anchor points during the iterative optimization normal equation construction stage; and reconstructing the increment of the normal equation using the second-order Laplace penalty term, so that the geometric reference anchor points act as the absolute geometric reference of the spatial scale during the adjustment process.

[0009] Preferably, the conformal mapping correction performed in step 5 includes the following sub-steps: Step 5-1, extracting the image reprojection residual distribution features corresponding to the low-weight image bundles within the connecting region; Step 5-2, constructing a local coordinate transformation operator using the object-space coordinates of the geometric reference anchor point; Step 5-3, using the local coordinate transformation operator to perform overall conformal translation, rotation, and scaling on the low-weight image bundles, so that the low-weight image bundles can achieve geometric stitching with the stable geometric skeleton in the object space while maintaining the internal relative topological structure.

[0010] Preferably, after step 5-3, the method further includes: real-time monitoring of the condition number of the equation matrix of the global bundle adjustment method; when the condition number tends to stabilize due to the intervention of the numerical damping factor, it is determined that the numerical oscillation of the connection region is suppressed.

[0011] Preferably, after step 4, the method further includes: monitoring the spatial distribution characteristics of the intersection constraint weights and reprojection residuals; extracting the structured residuals corresponding to the geometric reference anchor points; and adjusting the distortion correction coefficients in the interior orientation elements based on the structured residuals to compensate for the system errors introduced by model parameter drift.

[0012] Preferably, the identification of the connecting region also incorporates the image overlap distribution characteristics. When the image overlap is less than 50% and the weight gradient is greater than a preset threshold, the locking logic of the stable geometric skeleton is automatically activated.

[0013] Preferably, the process of generating intersection constraint weights also includes a compensation step for atmospheric refraction factor, which uses the photographic distance to linearly correct the optical path offset caused by atmospheric refraction, so as to improve the characterization accuracy of intersection constraint weights under long-distance photographic conditions.

[0014] Preferably, after step 5, the method further includes: outputting the real-world 3D model data after geometric closure processing, and smoothing the geometric patch connections of the real-world 3D model data in the connecting area to correct topological breaks.

[0015] Compared with the prior art, the beneficial effects of the present invention are: 1. In the photogrammetric space of multi-source image fusion, by establishing a deterministic mapping relationship between the image bundle stiffness coefficient and the intrinsic physical parameters of the sensor, the quantitative characterization of the spatial projection ray bundle constraint strength can be realized, which can effectively improve the condition number of the global adjustment method equation matrix, suppress the numerical oscillation caused by the difference in the physical properties of heterogeneous image sensors and the magnitude of shooting conditions, and ensure the logical self-consistency of observations at different scales under a unified weight framework.

[0016] 2. By utilizing the stiffness gradient changes between adjacent images to identify heterogeneous connection boundaries, and combining the local topological features extracted from high-stiffness images to construct a geometric skeleton, the system applies conformal constraints to low-stiffness rays within the overlapping region to lock the offset degrees of freedom of low-stiffness rays relative to the geometric skeleton. This eliminates geometric steps and topological breaks at the fusion interface of the cross-scale model, maintaining seamless stitching of the reconstructed model at the physical level. By monitoring the spatial orthogonal distribution characteristics of the image bundle stiffness coefficient and the reprojection residual, the system is able to distinguish between random observation noise and sensor model drift. When the structured residual corresponding to the high-stiffness image bundle is detected, the intrinsic parameter compensation vector is inverted using the residual gradient direction to achieve online compensation for sensor temperature drift under long-endurance operation, ensuring the continuous stability of reconstruction accuracy.

[0017] 3. By combining image bundle physical stiffness evaluation with dynamic damping adaptive adjustment, global drift is suppressed through the anchoring effect of high-stiffness local areas without relying on large-scale ground control points. This improves the convergence robustness of the adjustment system, ensures topological closure in the reconstruction model generation process under complex industrial conditions, and effectively suppresses nonlinear bending and model warping in long-range photogrammetric reconstruction. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the spatial reconstruction method of intersection constraint weights and stable geometric skeleton of the present invention; Figure 2 This invention provides a system architecture and data flow diagram that integrates multi-source data links and a high-performance reconfiguration engine. Detailed Implementation

[0019] The following disclosure is intended to explain and illustrate the present invention, and not to limit the scope of protection of the present invention. Unless otherwise specified, the technical features of the various embodiments and examples described below can be combined with each other. Those skilled in the art should understand that any modifications or equivalent substitutions made to the technical solutions of the present invention based on the technical essence determined by the present invention shall fall within the scope of protection of the present invention.

[0020] This invention provides a photogrammetric spatial reconstruction method based on multi-source image fusion. By establishing an intersection constraint strength model of image bundles in object space, it identifies the connection regions of heterogeneous images and extracts geometric reference anchor points. A stable geometric skeleton is constructed using differential regularization terms, and conformal mapping correction is performed on low-weight image bundles constrained by numerical damping factors, achieving geometric closure of multi-source images within a unified framework. After acquiring UAV-borne images and ground close-range images covering the target area, the system extracts the interior orientation elements and initial exterior orientation elements of each image acquisition unit. The interior orientation elements include focal length. In addition to internal geometric parameters such as pixel size δ, the initial exterior orientation elements are used to determine the position of the imaging center in the spatial coordinate system and the direction of the image optical axis. To quantitatively characterize the geometric constraint capability of the projection beams generated by different sensors on spatial points, the system analyzes the focal length f in the interior orientation elements and the imaging distance D and optical axis tilt angle θ corresponding to the initial exterior orientation elements, and calculates the intersection constraint weight w of each image beam in the spatial intersection topology according to the following formula: w=(f⋅cosθ) / (1000⋅D⋅δ), where w is the intersection constraint weight; f is the focal length in mm; θ is the optical axis tilt angle in ∘; D is the imaging distance in m; and δ is the pixel size in mm. During this calculation process, the system simultaneously uses the imaging distance D in combination with the standard atmospheric refraction model to perform equivalent correction on the projection optical path. Because long-distance aerial photography is affected by atmospheric density gradients, optical path deviations weaken the geometric constraint effectiveness of the image bundle. Through this physical environment compensation, it is ensured that the weight w can truly represent the geometric stiffness of image bundles of different scales in a unified object space, providing a logically consistent weighting benchmark for subsequent cross-scale adjustment. The intersection constraint weight w reflects the geometric stiffness of the projected ray bundle in the object space. In the initialization stage of the stochastic model of global bundle adjustment, the system uses the magnitude of the intersection constraint weight w as the initial input parameter of the observation weight matrix.

[0021] In overlapping regions of multi-source images, due to differences in the magnitude of sensor physical properties, the intersection constraint weights *w* between adjacent image bundles may fluctuate. The system monitors the spatial rate of change of these weights to determine the weight gradient ∇w. Specifically, a topological connectivity map of adjacent image bundles within the connecting region is established, and the first-order difference of the intersection constraint weights *w* between nodes in the topological connectivity map is calculated to obtain the weight gradient ∇w. When this gradient exceeds a preset topological mismatch threshold, the region is deemed to have a risk of numerical oscillation in the normal equations. The system extracts image bundles within the connecting region whose intersection constraint weights *w* are higher than the preset weight threshold as geometric reference anchor points. The selection of this preset threshold follows the statistical distribution characteristics of the weights across the entire field. By identifying the high-end statistical intervals of the weight distribution, the system dynamically determines a set of high-stiffness image bundles capable of supporting spatial scale benchmarks. This adaptive extraction mechanism avoids the failure of fixed thresholds during switching between different sensor sources, ensuring that the geometric skeleton constructed at the multi-source image interface has sufficient resistance to physical distortion and suppressing the risk of numerical oscillation caused by weight mismatch. The topological mismatch threshold T... ∇w The procedure adopts a statistical distribution calibration method. Fifty sets of heterogeneous image sequences covering the same overlapping area are collected in a static environment or standard calibration field. The distribution characteristics of the weight gradient ∇w for each image set with respect to the initial adjustment state are calculated. The root mean square error σ of the statistical weight gradient ∇w in the stationary observation area is determined. ∇w The topology mismatch threshold T ∇w Set to 3⋅σ ∇w The adjustment calculation monitors the geometric displacement deviation of the connection area and dynamically fine-tunes this value. When the geometric displacement deviation is greater than 1mm, it is adjusted by 0.5⋅σ. ∇w The threshold T is the step size decrease. ∇w The process continues until the geometric displacement deviation in the connecting area converges to within 0.5 mm. During the iterative optimization process of global bundle adjustment, the system applies a second-order differential regularization term based on the spatial position offset to the exterior orientation elements corresponding to the geometric reference anchor points. This procedure includes extracting the collinearity constraint equations formed by the geometric reference anchor points in the object space, and adding a second-order Laplace penalty term to the column vector of the exterior orientation elements corresponding to the geometric reference anchor points during the normal equation construction stage of iterative optimization. The second-order Laplace penalty term is used to reconstruct the increment of the normal equations, limiting the correction step size of the geometric reference anchor points in the iteration space, thereby constructing a stable geometric skeleton composed of geometric reference anchor points in the object space, so that it can act as the absolute geometric reference of the spatial scale during the adjustment process.

[0022] The system utilizes the intersection constraint weights *w* to construct a numerical damping factor, which is used to suppress numerical oscillations generated by low-weight image bundles during iteration. The physical meaning of the numerical damping factor is the correction weight coefficient of the iteration step size, and its value is positively correlated with the corresponding intersection constraint weights *w*. During the normal equation iterative update phase, the system establishes the diagonal element values ​​of the numerical damping factor by establishing a nonlinear mapping function between the intersection constraint weights *w* and the iteration step size. This function model sets the numerical damping factor as η = exp(-γ⋅(∇w / w)). Where η is the numerical damping factor, γ is the preset attenuation control coefficient, ∇w is the weight gradient, and w is the intersection constraint weight; this procedure uses the exponential attenuation feature to force the low-weight image bundle in the high gradient region to enter the minimum step size search state, thereby cutting off the error transmission path from the heterogeneous sensor noise to the stable geometric skeleton at the normal equation level; the regularization coefficient λ and the numerical damping factor η are co-tuned to determine the procedure, the normal equation is iteratively updated to monitor the condition number κ and eigenvalue distribution characteristics of the normal equation matrix in real time, and the regularization coefficient λ searches the logarithmic search interval

[10] . -2 10 1 Using a step size of 10 -1 Numerical scan, selecting a value that stabilizes the order of magnitude of the condition number κ by 10. 4 The following values, where the rate of change of the reprojection residual's descent slope is less than 0.05, are selected as preferred parameters. A numerical damping factor η = exp(-γ⋅(∇w / w)) is used to limit the step size of the low-weight image bundle within the iteration space. The attenuation control coefficient γ is chosen to minimize the difference in magnitude between the diagonal and off-diagonal elements of the normal equation matrix, cutting off the numerical oscillation path induced by differences in heterogeneous sensors and eliminating geometric steps at the image fusion interface. In the actual iteration process, the adjustment objective of the attenuation control coefficient γ is to optimize the condition number of the normal equation matrix. When the weight gradient ∇w between adjacent images is detected as potentially inducing matrix singularity, the damping factor cuts off the noise propagation path through nonlinear mapping. The system continuously monitors the eigenvalue distribution of the normal equation matrix until the matrix condition number stabilizes and the rate of change of the reprojection residual's descent slope meets the preset stationarity criterion. This ensures the geometric continuity of images at different resolutions at the fusion interface at the algorithm level.

[0023] The system constructs a local coordinate transformation operator with a stable geometric skeleton as the reference frame, forcing the low-weight image bundle to perform conformal mapping correction under the constraints of the local coordinate transformation operator. This is achieved by solving a seven-parameter similarity transformation model containing three rotational degrees of freedom, three translational degrees of freedom, and one scaling factor. In each adjustment iteration step, the reprojection residual center vector of all object points in the low-weight image bundle is extracted. Using the object coordinate true value of the geometric reference anchor point as the reference centroid, the rotation matrix is ​​solved by applying orthogonal matrix constraints, so that the low-weight image bundle as a whole achieves rigid alignment with the stable geometric skeleton without changing the relative angles of its internal projected ray bundles, eliminating the geometric step at the image fusion interface of the reconstructed model. After performing conformal mapping correction, the system monitors the condition number of the global bundle adjustment equation matrix in real time. When the condition number tends to stabilize, the numerical oscillation in the connection area is suppressed, and the system outputs geometrically closed real-scene 3D model data. The system smoothly connects geometric patches in the region, correcting topological breaks. To address sensor intrinsic parameter drift caused by environmental temperature differences, the system monitors the spatial distribution characteristics of the intersection constraint weights w and reprojection residuals, extracts the structured residuals corresponding to the geometric reference anchor points, and adjusts the distortion correction coefficients in the internal orientation elements based on the structured residual feedback when a structured distribution of reprojection residuals appears at the corresponding location of the high-stiffness image bundle, compensating for the system error introduced by model parameter drift. The core of this adjustment mechanism is to use the physical stability of the geometric reference anchor points as a locator. When the residuals of these high-stiffness image bundles show obvious radial or tangential vector patterns (i.e., structured residuals), the system determines that lens thermal deformation or parameter drift induced by environmental temperature differences has occurred. By inverting this structured component into an intrinsic parameter compensation vector, online compensation for changes in the physical state of the sensor is achieved, thereby suppressing model bending and warping in long-range reconstruction without relying on large-scale ground control points.

[0024] Example 1: In the application scenario of high-precision real-scene 3D reconstruction of a large cross-river bridge, photogrammetry operations face the objective constraint of cross-scale image fusion. The global frame image acquired through UAV oblique photography has a photographic distance D of 150m, while the detailed image acquired through a handheld device on the bridge pier has a photographic distance D shortened to 2m. This results in a mismatch in the geometric constraint strength of the spatial ray beams generated by different image sources. This asymmetry in physical properties leads to the degradation of the condition number of the normal equation matrix of the adjustment system, causing numerical oscillations at the model fusion interface. This manifests as a geometric step error of approximately 8cm at the scale transition point of the reconstruction model. The focal length f of the system acquisition device is 35mm, and the pixel size δ is 0.0039mm. After analyzing the interior orientation elements and initial exterior orientation elements of the image, the intersection constraint weights of the global image beams are calculated according to the formula w=(f⋅cosθ) / (1000⋅D⋅δ). Where w is the intersection constraint weight, f is the focal length in mm, θ is the optical axis tilt angle in ∘, D is the shooting distance in m, and δ is the pixel size in mm; the calculation results show that the intersection constraint weight w of the UAV high-altitude image bundle is around 60, while the intersection constraint weight w of the local close-up image bundle of the bridge pier reaches over 4400. By calculating the first difference of the intersection constraint weight w between adjacent image bundles in the connecting area, the weight gradient ∇w is determined, and the topological boundary where the weight magnitude changes drastically is identified, and the high-weight image bundles in it are locked as geometric reference anchor points.

[0025] During the iterative optimization process of global bundle adjustment, the system takes the collinearity constraint equations corresponding to the geometric reference anchor points as input. By adding a second-order Laplacian penalty term to the column vector of the exterior orientation elements corresponding to the geometric reference anchor points, the system achieves second-order differential regularization constraint on this part of the correction, thereby constructing a stable geometric skeleton in the object space composed of the geometric reference anchor points. The system uses the object space coordinates of the geometric reference anchor points to solve the local coordinate transformation operator, and forces the low-weight image bundle to perform overall conformal translation and rotation under the guidance of the local coordinate transformation operator. This completes the conformal mapping correction of the UAV image bundle at the fusion interface, and the processed reprojection residual distribution converges in the entire domain.

[0026] Example 2: Performed in a controlled industrial measurement laboratory environment, the stability of the reconstruction method was verified by constructing a cross-scale photogrammetry test platform, which includes a high-precision moving guide rail, a telephoto imaging unit, a close-range imaging unit, and a laser tracker with a positioning accuracy better than 0.015 mm. The telephoto imaging unit and the close-range imaging unit jointly observe the simulated bridge structure synchronously. The object coordinates of 12 geometric reference points on the surface of the simulated structure are measured by the laser tracker as a reference data stream for evaluating the reconstruction accuracy. The process of acquiring test data is achieved by controlling the position of the telephoto imaging unit on the guide rail to simulate the gradient change of the photographic distance D in the range of 20 m to 250 m.

[0027] During the parameter calibration phase, the sampling period is set based on a trade-off between the search complexity of image feature matching and physical motion displacement. When the displacement generated by the motion angular velocity of the imaging unit on the image plane exceeds 0.5 pixels, the sampling period is set to 200ms to suppress dynamic blur. To simulate atmospheric disturbances in a real photographic environment, Gaussian noise with a root mean square error of 0.3 pixels is superimposed on the initial observation values. The system analyzes the focal length f and the pixel size δ, and obtains the photographic distance D and the optical axis tilt angle θ through the initial adjustment values. According to the formula w=(f⋅cosθ) / (1000⋅D⋅δ), the intersection constraint weight w is calculated, where w is the intersection constraint weight, f is the focal length in mm, θ is the optical axis tilt angle in ∘, D is the photographic distance in m, and δ is the pixel size in mm. The system monitors the rate of change of the intersection constraint weight w between adjacent image bundles to determine the weight gradient ∇w, as shown in Table 1.

[0028] Table 1: Reconstruction accuracy data of the experimental group and the control group at different scale spans Control group A 8.0 512.4 1.00, meaning no damping applied. 4.82 Experimental group A 8.0 512.4 0.12 0.85 Control group B 40.0 2560.8 1.00, meaning no damping applied. 35.12 Experimental group B 40.0 2560.8 0.03 1.42 Control group C 100.0 6400.2 1.00, meaning no damping applied. 112.55 Experimental group C 100.0 6400.2 0.01 2.15 Boundary group 150.0 9600.5 0.005 6.84

[0029] Referring to the data in Table 1, with the ratio of photographic distance D ratio As the weight gradient ∇w increases, the control group shows a linear upward trend in D. ratio When the value reaches 100.0, due to the degradation of the condition number of the normal equation matrix, the geometric displacement deviation of the joint region reaches 112.55 mm, resulting in obvious topological breakage. The experimental group identified the joint interface by calculating the weight gradient ∇w and used the reciprocal of the intersection constraint weight w as a numerical damping factor η to limit the correction step size of the low-weight image bundle in the adjustment iteration, so that the geometric displacement deviation is always kept at a low level. ratio When the threshold value exceeds 150.0, the reconstruction bias shows a clear inflection point, indicating that the physical constraints provided by long-range images have approached the lower limit of the signal-to-noise ratio.

[0030] Example 3: This example combines Figures 1 to 2 This section explains the photogrammetric spatial reconstruction method based on multi-source image fusion, as follows: Figure 1As shown, step 1 acquires UAV-borne images and ground close-up images, and extracts the interior orientation elements and initial exterior orientation elements of each image acquisition unit. Step 2 analyzes the focal length, shooting distance, and optical axis tilt, establishes a geometric constraint strength model, and generates intersection constraint weights for each image bundle in the spatial intersection topology. Step 3 monitors the spatial change rate of the intersection constraint weights to determine the weight gradient, identifies the connection region, and extracts high-weight image bundles as geometric reference anchor points. Based on this, step 4 applies a second-order differential regularization term based on spatial position offset to the geometric reference anchor points, limits the correction step size, and constructs a stable geometric skeleton in the object space. Finally, step 5 uses a numerical damping factor to suppress oscillations and forces low-weight image bundles to perform conformal mapping correction through local coordinate transformation to eliminate the geometric step at the fusion interface.

[0031] like Figure 2 As shown, the system architecture includes a data acquisition terminal layer, which encompasses an aerial drone platform equipped with a high-resolution camera, a ground measurement system including a ground close-range camera and a mobile terminal, and an external image data source that accesses satellite remote sensing images. These devices are connected to an image data center server via a multi-source image data transmission link. The server is equipped with a large-capacity storage array and internally deploys a raw multi-source image database as well as an intermediate process and parameter database for storing weight maps and skeleton data. It is connected to a high-performance computing server or cluster via a high-speed data I / O channel for raw image reading and intermediate result access. This computing unit uses GPU acceleration and multi-core parallel processing technology to run a photogrammetric spatial reconstruction engine. Internally, it includes a beam constraint strength modeling module, a heterogeneous region identification and classification module, a stable geometric skeleton modeling module using differential regularization, a conformal mapping correction module applying weighted damping constraints, and a geometric closure and seamless stitching module. Finally, the reconstruction results are output to a professional application workstation via a reconstruction result publishing link. This workstation serves as a high-performance graphics display terminal, equipped with 3D visualization and analysis software and a human-computer interaction interface.

[0032] Example 4: In the digital flexible assembly scenario of large aerospace composite panels, the photogrammetry system faces the failure of physical displacement constraints caused by environmental temperature fluctuations. When the ambient temperature in the assembly plant shifts by more than 8°C within 12 hours, the images acquired by the fixed observation station on the top of the plant and the local assembly position observation station show significant differences in the distance between the image centers and the physical pixel resolution. This causes the intersection constraint weight w at the interface of the heterogeneous images to experience a gradient jump of more than 2000 orders of magnitude, inducing numerical singularities in the global bundle adjustment equation matrix. This results in a normal displacement deviation of the reconstructed model at the panel joint reaching the order of 2.15mm, exceeding the 0.5mm allowable range for engineering assembly. The system executes a parameter calibration procedure based on statistical characteristics to determine the topology mismatch threshold T. ∇wThe quantization boundary values ​​include acquiring 50 consecutive heterogeneous image sequences covering the same overlapping area under static conditions, and calculating the weight gradient of each image pair under the initial adjustment state. Distribution characteristics were analyzed, and the root mean square error σ of the weight gradient ∇w in the stationary observation region was obtained using computational software. ∇w And set the topology mismatch threshold T ∇w Set to 3⋅σ ∇w Real-time monitoring of the weight gradient ∇w between adjacent image bundles; when the observed weight gradient ∇w reaches the calibrated topological mismatch threshold T... ∇w When the time is right, the system determines that the user has entered the heterogeneous image fusion interface and automatically activates the geometric skeleton locking operator.

[0033] In the iterative optimization procedure for global bundle adjustment, the system embeds a discretized second-order differential penalty term into the coefficient matrix of the normal equations. This is achieved by denoting the column vector of the exterior orientation elements of the geometric reference anchor points as Δx, and accumulating the second-order differential operator in the corresponding diagonal matrix unit. This operator follows the mathematical constraint relation (B... T PB+λL)Δx=B T Pl, where B is the partial derivative matrix of the error equation, P is the observation weight matrix, and Δx is the increment of the exterior orientation element to be determined. Let be the column vector of observation residuals, λ be the regularization coefficient, and L be the second-order Laplacian matrix constructed from the topological connections of the geometric reference anchor points. This matrix L is used to restrict the non-physical displacements of the geometric reference anchor points in the iteration space, enabling the high-stiffness ray beam to form a stable geometric skeleton with physical distortion resistance in the object space. To address sensor intrinsic parameter drift caused by environmental temperature fluctuations, the system executes a compensation program based on the distribution characteristics of the reprojection residual vector, extracts the reprojection residual column vector r corresponding to the geometric reference anchor points, and calculates the spatial cross-correlation coefficient between the radial component of the reprojection residual column vector r and the corresponding intersection constraint weight w. When the gradient vector of the reprojection residual is greater than 0.85 and points to the principal point of the image, it is determined that sensor model parameter drift caused by lens thermal deformation has occurred. The system extracts the mean offset of the structured residual and maps it to an intrinsic parameter correction vector Δk. By introducing the intrinsic parameter correction vector Δk into the distortion correction model of the collinear equation, the principal point offset and lens distortion term changes are offset, thus realizing online compensation for changes in the physical environment. Assembly monitoring data shows that after the system executes the above processing procedure, the panel normal reconstruction error is reduced from 2.15mm to 0.18mm, and the condition number of the adjustment system normal equation matrix is ​​reduced from 10. 8 Reduced to 10 3 This eliminates numerical oscillations caused by the mismatch of sensor physical properties, maintains the topological continuity of geometric patches in the real-scene 3D model at the image fusion interface, and corrects the geometric steps generated in cross-scale spatial reconstruction.

[0034] Example 5: In the on-site deployment process of the real-scene reconstruction project of a super-large water conservancy hub, the system starts the pre-calibration program to correct the deviation between the nominal parameters of the sensor and the physical conditions. The program runs in an environment with relative humidity below 65% and wind speed less than 3m / s. By setting up 20 sets of calibration points with known object coordinates covering the edge of the reconstruction area, the system acquires calibration images of each imaging unit at gradient height and different angular rates, calculates the effective focal length f of the sensor and the geometric correction coefficient of the pixel size δ under the current temperature field, and establishes the Euclidean distance baseline between the initial position point of the station and the center point of the calibration field to calibrate the measurement zero point of the initial photographic distance D. The parameter matrix determined in this way serves as the reference system for generating the intersection constraint weights w.

[0035] Before starting formal observation data acquisition, the system executes an offline injection program for the weight gradient benchmark, statistically analyzes the distribution characteristics of the weight gradient ∇w of different imaging units under standard overlap, calculates the intersection angle components of each image bundle in the overlapping area to determine the initial weight distribution matrix, and introduces a regularization term with physical constraint directionality in the initial value calculation stage of the adjustment method equation to suppress spurious weight oscillations caused by the coupling of atmospheric refraction and sensor thermal noise, so that the condition number of the global normal equation matrix reaches a stable state within the first 5 iterations after system startup, and the scale consistency deviation of the reconstructed object points at the fusion boundary is on the order of millimeters.

[0036] Example 6: During the on-site deployment of the safety inspection of high-speed railway bridge pier structures, the system executes the calibration procedure of the sensor stiffness reference matrix. It is started under conditions where the relative humidity is below 65% and the equipment vibration frequency caused by wind speed avoids the 50Hz resonant frequency. By deploying 15 sets of corresponding control points with known object coordinates within a local area, images of each imaging unit at gradient heights are acquired, and the geometric correction coefficients of the effective focal length f and pixel size δ of the sensor under the current temperature field are calculated. Gradient imaging is performed within the shooting distance D from 1.5m to 30m to extract the distribution reference of the intersection constraint weights w. The eigenvalue distribution of the normal equation matrix is ​​monitored in real time, and the matrix condition number κ is calculated. The system is then executed within the logarithmic search interval of the regularization coefficient λ, with a step size of 10. -1 Numerical scans were performed to determine the order of magnitude of the matrix condition number κ to be stabilized at 10. 4 The regularization coefficient λ is defined as follows, where the rate of change of the slope of the reprojection residual is less than 0.05.

[0037] During the refactoring process, the system utilizes the calibrated regularization coefficients. The online adaptive adjustment procedure is initiated, and the weight gradient ∇w data between adjacent images are collected in real time. When the displacement of the imaging unit caused by environmental disturbance is determined by the intrinsic parameter drift discrimination logic to produce unstructured residuals, the second-order Laplacian matrix constructed from the topological connection relationship of the geometric reference anchor points is extracted. Furthermore, the local coordinate transformation operator is recalculated within each step of the adjustment iteration, forcing the low-weight image bundle to perform conformal mapping correction under the physical constraints of a stable geometric skeleton. Statistical data shows that, under the condition of a cross-scale ratio of 20, the geometric displacement deviation of the connecting area is reduced from 3.85mm to 0.42mm, and the reconstructed object points maintain the topological continuity of the physical scale at the image fusion interface.

[0038] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the present invention can be implemented in other specific forms without departing from the spirit or essential characteristics of the present invention.

[0039] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention.

Claims

1. A photogrammetric spatial reconstruction method based on multi-source image fusion, characterized in that, Includes the following steps: Step 1: Acquire UAV-borne images and close-up ground images, and extract the interior orientation elements and initial exterior orientation elements of each image acquisition unit; Step 2: Analyze the focal length in the interior orientation elements and the photographic distance and optical axis tilt angle corresponding to the initial exterior orientation elements. Based on the proportional relationship between focal length and photographic distance and combined with the cosine component of optical axis tilt angle, establish a geometric constraint strength model that characterizes the intersection stability of image bundles in object space, so as to generate the intersection constraint weights of each image bundle in the spatial intersection topology. Step 3: Monitor the spatial change rate of the intersection constraint weights between adjacent image bundles to determine the weight gradient. Based on the weight gradient, identify the connection area between the UAV-borne image and the ground close-range image, and extract the image bundles in the connection area whose intersection constraint weights are higher than a preset threshold as geometric reference anchor points. Step 4: In the iterative optimization process of global bundle adjustment, a second-order differential regularization term based on the spatial position offset is applied to the exterior orientation element corresponding to the geometric reference anchor point. The second-order differential regularization term is used to limit the correction step size of the geometric reference anchor point in the iteration space, thereby constructing a stable geometric skeleton composed of geometric reference anchor points in the object space. Step 5: Construct a numerical damping factor using intersection constraint weights. Use the numerical damping factor to suppress the numerical oscillations generated by the low-weight image bundle during the iteration process. By constructing a local coordinate transformation operator with a stable geometric skeleton as the reference system, force the low-weight image bundle to perform conformal mapping correction under the constraint of the local coordinate transformation operator, thereby locking the geometric offset degree of freedom of the low-weight image bundle relative to the stable geometric skeleton, so as to eliminate the geometric step of the reconstruction model at the image fusion interface.

2. The photogrammetric spatial reconstruction method based on multi-source image fusion according to claim 1, characterized in that, The logical rule followed in step 2 for calculating the intersection constraint weights w of the image bundle is as follows: obtain the focal length f, the shooting distance D, the optical axis tilt angle θ, and the pixel size δ of the image acquisition unit, and calculate the intersection constraint weights w according to the following formula: w=(f⋅cosθ) / (1000⋅D⋅δ), where w is the intersection constraint weight, f is the focal length in mm, θ is the optical axis tilt angle in ∘, D is the shooting distance in m, and δ is the pixel size in mm; and in the initialization stage of the stochastic model of global bundle adjustment, the magnitude of the normalized intersection constraint weights w is used as the initial input parameter of the observation weight matrix.

3. The photogrammetric spatial reconstruction method based on multi-source image fusion according to claim 1, characterized in that, Step 3 includes the following sub-steps: Step 3-1, using the existence of shared feature points between image bundles as the basis for adjacency determination, establish a topological connection map between adjacent image bundles in the connection area; Step 3-2, calculate the first-order difference of the intersection constraint weights between each node in the topological connection map to obtain the weight gradient; Step 3-3, when the weight gradient exceeds the preset topological mismatch threshold, determine that there is a risk of numerical oscillation of the normal equation induced by the difference in the magnitude of the physical properties of the sensors in the connection area.

4. The photogrammetric spatial reconstruction method based on multi-source image fusion according to claim 1, characterized in that, Step 4 involves applying a second-order differential regularization term based on spatial position offset, which includes: extracting the collinearity constraint equations formed by the geometric datum anchor points in the object space; adding a second-order Laplace penalty term to the column vector of the exterior orientation elements corresponding to the geometric datum anchor points during the iterative optimization of the normal equation construction stage; and reconstructing the increment of the normal equations using the second-order Laplace penalty term, so that the geometric datum anchor points act as the absolute geometric datum of the spatial scale during the adjustment process.

5. The photogrammetric spatial reconstruction method based on multi-source image fusion according to claim 1, characterized in that, Step 5 includes the following sub-steps for performing conformal mapping correction: Step 5-1, extract the image reprojection residual distribution features corresponding to the low-weight image bundles in the connecting region; Step 5-2, construct a local coordinate transformation operator using the object-space coordinates of the geometric reference anchor point; Step 5-3, use the local coordinate transformation operator to perform overall conformal translation, rotation and scaling on the low-weight image bundles.

6. The photogrammetric spatial reconstruction method based on multi-source image fusion according to claim 5, characterized in that, Step 5-3 is followed by: real-time monitoring of the condition number of the equation matrix of the global bundle adjustment method; when the condition number tends to stabilize due to the intervention of the numerical damping factor, it is determined that the numerical oscillation of the connection region is suppressed.

7. The photogrammetric spatial reconstruction method based on multi-source image fusion according to claim 1, characterized in that, Step 4 is followed by: monitoring the spatial distribution characteristics of intersection constraint weights and reprojection residuals; extracting the structured residuals corresponding to the geometric reference anchor points; and adjusting the distortion correction coefficients in the interior orientation elements based on the structured residuals to compensate for the system errors introduced by model parameter drift.

8. The photogrammetric spatial reconstruction method based on multi-source image fusion according to claim 1, characterized in that, The identification of the connecting region also incorporates the characteristics of image overlap distribution. When the image overlap is less than 50% and the weight gradient is greater than the preset threshold, the locking logic of the stable geometric skeleton is automatically activated.

9. The photogrammetric spatial reconstruction method based on multi-source image fusion according to claim 1, characterized in that, The process of generating intersection constraint weights also includes a compensation step for atmospheric refraction factor. The optical path offset caused by atmospheric refraction is linearly corrected by using the photographic distance to improve the representation accuracy of intersection constraint weights under long-distance photographic conditions.

10. A photogrammetric spatial reconstruction method based on multi-source image fusion according to claim 1, characterized in that, Step 5 is followed by: outputting the real-world 3D model data after geometric closure processing, and smoothing the geometric patch connections of the real-world 3D model data in the connecting areas to correct topological breaks.

Images