Land space image analysis method based on photogrammetry

Abstract

The present application relates to the field of photogrammetry, and discloses a land space image analysis method based on photogrammetry, comprising: acquiring aerial images and attitude data of a survey area to generate sparse point clouds; distinguishing rigid and non-rigid feature points based on semantic segmentation results; establishing a local vertical reference system for non-rigid points and configuring anisotropic parameters with transverse variance greater than longitudinal variance; mapping the parameters to a global weight matrix using a coordinate rotation matrix, and substituting them into a bundle adjustment model to complete iterative solution, the present application effectively realizes the automatic suppression of non-rigid disturbance of vegetation and water area by constructing an anisotropic weight determination mechanism coupled with physical properties and geometric constraints, and improves the geometric accuracy and stability of three-dimensional reconstruction in complex land space scenes.

Description

Technical Field

[0001] This invention belongs to the field of photogrammetry technology, and in particular relates to a method for analyzing land and space images based on photogrammetry. Background Technology

[0002] In current aerial and aerospace photogrammetry engineering practice, bundle adjustment is the core solution method for constructing high-precision 3D geographic information models. This method is based on the collinearity equation theory and solves the image exterior orientation elements and the 3D coordinates of ground feature points by minimizing the reprojection error. Under ideal rigid geometric scenarios, the least squares adjustment algorithm based on the Gauss-Markov model can achieve pixel-level measurement accuracy and is widely used in topographic mapping, urban 3D modeling and land monitoring.

[0003] Current technologies for improving the geometric accuracy of measurement results often focus on enhancing the mechanical pointing accuracy and automation level of hardware systems to avoid errors. Few consider the constraints imposed on adjustment models by the non-rigid physical characteristics of ground features at the algorithmic level. For example, Chinese invention patent CN103837138B discloses a precision photogrammetry robot that constructs an automated system integrating remote sensing, telemetry, and 3D attitude control. Through high-precision motor drives and precision ranging units, it achieves automatic target aiming and 3D coordinate telemetry. While this technology, with its hardware upgrades, has advantages in fixed-point measurement of rigid targets, its technical concept is still based on the assumption that the observed object is a rigid, static body. It lacks the ability to statistically model dynamic changes in the observation environment. When applied to land space analysis involving large amounts of vegetation or water, this purely hardware-accuracy-dependent approach cannot distinguish between sensor noise and anisotropic displacements of ground features caused by wind loads. Furthermore, when processing non-rigid feature points, it cannot adaptively relax the level through mathematical models. Constraints make it difficult to meet the accuracy requirements of high-fidelity 3D reconstruction in complex scenarios. When conducting detailed analysis of land space containing large areas of vegetation cover, complex water areas, and high-density urban building clusters, there is a contradiction between the complexity of the physical environment and the idealized assumptions of traditional mathematical models. Traditional bundle adjustment is usually based on the assumptions of homogeneous rigid bodies and isotropic errors, assuming that the feature points in the observation scene are all geometrically fixed rigid points, and that their observation errors follow the same statistical distribution in all directions of 3D space. In actual engineering, this static and isotropic mathematical description is difficult to adapt to the real surface environment that presents physical rheological characteristics. Constrained by the gravitational field and biomechanical structure, vegetation feature points usually have high positional stability in the vertical direction, but in the horizontal direction, they are affected by wind loads and produce large-scale random swaying. The mirror reflection virtual image formed by the glass curtain wall of high-rise buildings presents rigid features in texture, but does not meet the same source intersection condition in the multi-view geometric light path.

[0004] Therefore, the technical problem to be solved by this invention is how to construct a bundle adjustment model that can accurately map the physical anisotropy characteristics of ground objects and dynamically identify the geometric and physical reliability of observation data, so as to achieve high-fidelity three-dimensional spatial analysis under complex land cover conditions. Summary of the Invention

[0005] This invention provides a method for analyzing land spatial images based on photogrammetry, comprising the following steps: The system acquires aerial photographic image sequences of the area to be tested, as well as spatiotemporally synchronized positioning and attitude data with the aerial photographic image sequences. It then performs feature extraction and image matching on the aerial photographic image sequences to generate an initial sparse point cloud containing three-dimensional spatial coordinates. Perform semantic segmentation of ground cover categories on aerial photographic image sequences, and divide the feature points in the initial sparse point cloud into a set of rigid feature points and a set of non-rigid feature points based on the segmentation results. For each feature point in the set of non-rigid feature points, establish a local perpendicular reference frame with the direction of the gravity vector as the Z-axis; Within a local vertical reference frame, a longitudinal variance component along the Z-axis and a transverse variance component along a plane perpendicular to the Z-axis are configured for the feature points, wherein the value of the transverse variance component is set to be greater than the value of the longitudinal variance component, thereby constructing an anisotropic variance diagonal matrix. The anisotropic variance diagonal matrix is ​​transformed from the local perpendicular reference system to the photogrammetric global coordinate system using a coordinate rotation matrix. The anisotropic covariance matrix in off-diagonal form is calculated, and the anisotropic weight matrix is ​​obtained by inverting the anisotropic covariance matrix. The anisotropic weight matrix is ​​substituted into the error equation of the bundle adjustment model as the stochastic model parameters of the observations. The least squares iterative adjustment is performed to output the adjusted exterior orientation elements and the three-dimensional coordinates of the ground points.

[0006] Preferably, the step of establishing a local perpendicular reference system with the gravity vector direction as the Z-axis specifically includes: extracting gravity acceleration vector data from the positioning and attitude data, or calculating the fitting plane normal vector of the initial sparse point cloud in the neighborhood of the current feature point, and defining the obtained vector direction as the gravity vector direction; using the gravity vector direction as the Z-axis of the local perpendicular reference system, and constructing two orthogonal vectors perpendicular to the Z-axis as the X-axis and Y-axis of the local perpendicular reference system, respectively; calculating the direction cosine matrix of the photogrammetric global coordinate system relative to the local perpendicular reference system, and determining the direction cosine matrix as the coordinate rotation matrix.

[0007] Preferably, the step of configuring a longitudinal variance component along the Z-axis and a lateral variance component along a plane perpendicular to the Z-axis for the feature points specifically includes: identifying the specific land cover category of the non-rigid feature points; if the land cover category is vegetation, then setting the longitudinal variance component to a first variance threshold and setting the lateral variance component to a second variance threshold, wherein the second variance threshold is at least 100 times the first variance threshold, so as to reduce the constraint weight on the horizontal displacement of vegetation feature points in the adjustment solution; if the land cover category is water, then setting the longitudinal variance component to the first variance threshold and maintaining the elevation constraint on water surface feature points.

[0008] Preferably, the calculation process of the anisotropic covariance matrix satisfies the following mathematical relationship: ,in, R is the anisotropic covariance matrix, and R is the coordinate rotation matrix. It is the transpose of the coordinate rotation matrix. Let P be the anisotropic variance diagonal matrix; the formula for calculating the anisotropic weight matrix P is: .

[0009] Preferably, the step of performing least squares iterative adjustment calculation further includes a weighting step based on the regional geometric rigidity index: calculating the regional geometric rigidity index of the initial sparse point cloud in the local region, the regional geometric rigidity index being used to characterize the ratio of the distribution density of rigid feature points to non-rigid feature points in the current adjustment unit; determining the dynamic balance coefficient based on the regional geometric rigidity index; adjusting the weight ratio of the image observation equation and the auxiliary navigation observation equation in the bundle adjustment model using the dynamic balance coefficient; and increasing the weight ratio of positioning and attitude determination data in the adjustment calculation when the regional geometric rigidity index is lower than a preset stability threshold using the dynamic balance coefficient.

[0010] Preferably, the calculation logic for the regional geometric stiffness index is as follows: a topological neighborhood is constructed centered on the currently solved image frame, and the number of rigid feature points within this neighborhood is counted. The number of non-rigid feature points ; Calculate the proportion of rigid points The calculation formula is: ; the proportion of rigid points As a regional geometric stiffness index, among which The value of is positively correlated with the regional geometric stiffness index.

[0011] Preferably, the step of performing least squares iterative adjustment solution specifically includes: calculating the reprojection residual vector of each feature point in a single iteration; calculating the equivalent weight factor using a robust estimation function based on the anisotropic weight matrix and the reprojection residual vector; iteratively updating the anisotropic weight matrix using the equivalent weight factor to reduce the influence of feature points whose reprojection residual vector exceeds three times the mean square error limit on the adjustment results, until the error equation converges.

[0012] Preferably, the step of performing semantic segmentation of land cover categories on the aerial photographic image sequence specifically includes: inputting the aerial photographic image sequence into a pre-trained semantic segmentation neural network to generate a pixel-level semantic mask; performing spatial projection matching between the pixel-level semantic mask and the initial sparse point cloud to assign a unique semantic label to each three-dimensional point in the initial sparse point cloud; performing point cloud filtering based on the semantic labels to remove feature points labeled as dynamically moving objects, retaining feature points labeled as buildings and roads as rigid feature points, and retaining feature points labeled as woodlands and water surfaces as non-rigid feature points.

[0013] Preferably, the method further includes: using the exterior orientation elements optimized by least squares iterative adjustment and the three-dimensional coordinates of ground points to perform dense matching processing on the aerial photographic image sequence to generate a high-density point cloud; constructing a digital surface model based on the high-density point cloud, and using the aerial photographic image sequence to perform texture mapping on the digital surface model to generate a real-scene three-dimensional model.

[0014] Compared with existing technologies, the land spatial image analysis method based on photogrammetry of this invention has the following advantages: 1. In the construction of the stochastic model for bundle adjustment, an anisotropic weight matrix strictly aligned with the local gravity vector is established. Taking advantage of the vertical stability of semi-rigid features such as vegetation due to the physical and mechanical constraints of gravity and growth structure, the variance components of its vertical degrees of freedom are specifically locked in the adjustment solution, while the variance constraints of its horizontal degrees of freedom affected by wind load are relaxed. This mathematical processing logic enables vegetation feature points, which are considered gross errors or noise in traditional measurements, to be transformed into effective elevation control constraints. Thus, without relying on external ground control points, the accurate transfer and convergence of elevation benchmarks for forest or vegetation-covered areas can be achieved by relying on the data redundancy of the image itself, avoiding the problems of voids or rank deficiencies in the elevation model caused by the forced removal of vegetation points.

[0015] 2. Based on semantic-assisted weighting, a geometric homology check mechanism based on the multi-view ray intersection discreteness is introduced. Utilizing the physical optical law that real rigid points satisfy the collinearity equation intersection condition, while specular reflection virtual images or dynamic interference points inevitably diverge in the object space, a second layer of objective verification defense independent of texture features is constructed. For glass curtain wall reflection points or water surface light spots that are misjudged as rigid objects by the semantic segmentation algorithm, once their ray intersection discreteness exceeds the geometric homology threshold, the system automatically triggers the weighted circuit breaking logic to exclude them from the solution mesh, blocking the propagation path of systematic model errors. This ensures that in urban canyon environments with strong reflections or repetitive textures, the 3D reconstruction results strictly converge to the real physical surface, eliminating the distortion and artifact phenomena commonly found in building models.

[0016] 3. By utilizing the geometric rigidity distribution characteristics of feature points across the entire map, a scene rigidity entropy index is constructed. This index is then used as an adjustment lever to dynamically modulate the weight ratio of exterior orientation element observations and visual feature observations in the adjustment equation. In urban areas where satellite positioning signals are susceptible to multipath interference but visual textures are rich, the weight of exterior orientation elements is automatically suppressed, and the rigidity of buildings is used to lock the network shape to resist positioning drift. In contrast, in areas with poor visual features, such as water bodies, the weight of inertial navigation data is automatically increased to prevent solution divergence. This endogenous adjustment mechanism, which does not require external filters, ensures that the system can always automatically anchor to observation sources with higher physical reliability during continuous operations with different land cover types, guaranteeing the spatial consistency of geometric accuracy of long-baseline, large-scale mapping results. Attached Figure Description

[0017] Figure 1 This is a flowchart of the image analysis process that integrates semantic segmentation and anisotropic weighting according to the present invention. Figure 2 This is the dynamic weighting logic diagram of the regional geometric stiffness index of this invention. Detailed Implementation

[0018] The technical solutions of the embodiments of this application will be clearly described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of this application are within the scope of protection of this application.

[0019] It should be noted that all directional and positional terms used in this invention, such as: up, down, left, right, front, back, vertical, horizontal, inner, outer, top, bottom, transverse, longitudinal, center, etc., are only used to explain the relative positional relationship and connection between components in a specific state (as shown in the accompanying drawings). They are only for the convenience of describing this invention and do not require that this invention be constructed and operated in a specific orientation. Therefore, they should not be construed as limiting this invention. In addition, the descriptions of "first," "second," etc., in this invention are for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated.

[0020] In the description of this invention, unless otherwise explicitly specified and limited, the terms installation, connection, and linking should be interpreted broadly. For example, they can refer to fixed connections, detachable connections, or integral connections; they can refer to mechanical connections; they can refer to direct connections or indirect connections through an intermediate medium; they can refer to the internal connection of two components. For those skilled in the art, the specific meaning of the above terms in this invention can be understood according to the specific circumstances.

[0021] In the description of this specification, references to the terms "an embodiment," "some embodiments," "illustrative embodiments," "examples," "specific examples," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the present invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example, and the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0022] A method for analyzing land spatial images based on photogrammetry includes the following steps: The system acquires aerial photographic image sequences of the area to be tested, as well as spatiotemporally synchronized positioning and attitude data with the aerial photographic image sequences. It then performs feature extraction and image matching on the aerial photographic image sequences to generate an initial sparse point cloud containing three-dimensional spatial coordinates. Perform semantic segmentation of ground cover categories on aerial photographic image sequences, and divide the feature points in the initial sparse point cloud into a set of rigid feature points and a set of non-rigid feature points based on the segmentation results. For each feature point in the set of non-rigid feature points, establish a local perpendicular reference frame with the direction of the gravity vector as the Z-axis; Within a local vertical reference frame, a longitudinal variance component along the Z-axis and a transverse variance component along a plane perpendicular to the Z-axis are configured for the feature points, wherein the value of the transverse variance component is set to be greater than the value of the longitudinal variance component, thereby constructing an anisotropic variance diagonal matrix. The anisotropic variance diagonal matrix is ​​transformed from the local perpendicular reference system to the photogrammetric global coordinate system using a coordinate rotation matrix. The anisotropic covariance matrix in off-diagonal form is calculated, and the anisotropic weight matrix is ​​obtained by inverting the anisotropic covariance matrix. The anisotropic weight matrix is ​​substituted into the error equation of the bundle adjustment model as the stochastic model parameters of the observations. The least squares iterative adjustment is performed to output the adjusted exterior orientation elements and the three-dimensional coordinates of the ground points.

[0023] Preferably, the step of establishing a local perpendicular reference system with the gravity vector direction as the Z-axis specifically includes: extracting gravity acceleration vector data from the positioning and attitude data, or calculating the fitting plane normal vector of the initial sparse point cloud in the neighborhood of the current feature point, and defining the obtained vector direction as the gravity vector direction; using the gravity vector direction as the Z-axis of the local perpendicular reference system, and constructing two orthogonal vectors perpendicular to the Z-axis as the X-axis and Y-axis of the local perpendicular reference system, respectively; calculating the direction cosine matrix of the photogrammetric global coordinate system relative to the local perpendicular reference system, and determining the direction cosine matrix as the coordinate rotation matrix.

[0024] Preferably, the step of configuring a longitudinal variance component along the Z-axis and a lateral variance component along a plane perpendicular to the Z-axis for the feature points specifically includes: identifying the specific land cover category of the non-rigid feature points; if the land cover category is vegetation, then setting the longitudinal variance component to a first variance threshold and setting the lateral variance component to a second variance threshold, wherein the second variance threshold is at least 100 times the first variance threshold, so as to reduce the constraint weight on the horizontal displacement of vegetation feature points in the adjustment solution; if the land cover category is water, then setting the longitudinal variance component to the first variance threshold and maintaining the elevation constraint on water surface feature points.

[0025] Preferably, the calculation process of the anisotropic covariance matrix satisfies the following mathematical relationship: ,in, R is the anisotropic covariance matrix, and R is the coordinate rotation matrix. It is the transpose of the coordinate rotation matrix. Let P be the anisotropic variance diagonal matrix; the formula for calculating the anisotropic weight matrix P is: .

[0026] Preferably, the step of performing least squares iterative adjustment calculation further includes a weighting step based on the regional geometric rigidity index: calculating the regional geometric rigidity index of the initial sparse point cloud in the local region, the regional geometric rigidity index being used to characterize the ratio of the distribution density of rigid feature points to non-rigid feature points in the current adjustment unit; determining the dynamic balance coefficient based on the regional geometric rigidity index; adjusting the weight ratio of the image observation equation and the auxiliary navigation observation equation in the bundle adjustment model using the dynamic balance coefficient; and increasing the weight ratio of positioning and attitude determination data in the adjustment calculation when the regional geometric rigidity index is lower than a preset stability threshold using the dynamic balance coefficient.

[0027] Preferably, the calculation logic for the regional geometric stiffness index is as follows: a topological neighborhood is constructed centered on the currently solved image frame, and the number of rigid feature points within this neighborhood is counted. The number of non-rigid feature points ; Calculate the proportion of rigid points The calculation formula is: ; the proportion of rigid points As a regional geometric stiffness index, among which The value of is positively correlated with the regional geometric stiffness index.

[0028] Preferably, the step of performing least squares iterative adjustment solution specifically includes: calculating the reprojection residual vector of each feature point in a single iteration; calculating the equivalent weight factor using a robust estimation function based on the anisotropic weight matrix and the reprojection residual vector; iteratively updating the anisotropic weight matrix using the equivalent weight factor to reduce the influence of feature points whose reprojection residual vector exceeds three times the mean square error limit on the adjustment results, until the error equation converges.

[0029] Preferably, the step of performing semantic segmentation of land cover categories on the aerial photographic image sequence specifically includes: inputting the aerial photographic image sequence into a pre-trained semantic segmentation neural network to generate a pixel-level semantic mask; performing spatial projection matching between the pixel-level semantic mask and the initial sparse point cloud to assign a unique semantic label to each three-dimensional point in the initial sparse point cloud; performing point cloud filtering based on the semantic labels to remove feature points labeled as dynamically moving objects, retaining feature points labeled as buildings and roads as rigid feature points, and retaining feature points labeled as woodlands and water surfaces as non-rigid feature points.

[0030] Preferably, the method further includes: using the exterior orientation elements optimized by least squares iterative adjustment and the three-dimensional coordinates of ground points to perform dense matching processing on the aerial photographic image sequence to generate a high-density point cloud; constructing a digital surface model based on the high-density point cloud, and using the aerial photographic image sequence to perform texture mapping on the digital surface model to generate a real-scene three-dimensional model.

[0031] Example 1: In an aerial photogrammetry operation targeting complex land areas with dense vegetation cover and water boundaries, the data processing system receives a multi-view overlapping image set acquired by an airborne five-lens oblique photogrammetry camera, along with its spatiotemporally synchronized positioning and attitude data. Under continuous crosswind conditions, the system performs feature extraction and image matching based on the SIFT operator on the image set to generate an initial sparse point cloud. A pixel-level semantic mask is then matched with the spatial projection of the initial sparse point cloud. Specifically, the semantic mask generation process involves using a ResNet-101 deep neural network that has undergone transfer learning on the UAVid UAV oblique photogrammetry dataset. The input image is uniformly cropped to a resolution of 1024 x 1024 pixels. The network output layer, after passing through a Softmax classifier, outputs a four-channel probability map containing vegetation, water bodies, buildings, and roads. During the projection matching stage, a depth buffer mechanism is introduced to eliminate semantic mapping deviations caused by field-of-view occlusion. The depth values ​​of multiple 3D feature points projected to the same pixel coordinates are calculated in the current camera coordinate system. The system sorts the data by depth and maps semantic mask labels only to feature points with the smallest depth values. Label transfer is performed when the difference between the reprojection depth of the feature point and the pixel value corresponding to the depth map generated by dense matching is less than the preset projection tolerance of 0.2 meters, so that the semantic category is associated with the visible feature points on the physical surface.

[0032] For non-rigid feature points in the initial sparse point cloud identified as vegetation categories by a semantic segmentation network, the system reads the instantaneous gravity acceleration vector data recorded in the positioning and attitude data, standardizes it, and defines it as the Z-axis direction of the local vertical reference frame at the feature point. Then, the system calls upon the instantaneous wind speed vector recorded by the airborne meteorological sensor, projects it onto a horizontal plane perpendicular to the Z-axis, defines the projection direction as the X-axis (i.e., downwind direction) of the local vertical reference frame, and defines the direction perpendicular to both the X-axis and Z-axis as the Y-axis (i.e., crosswind direction), thus constructing a right-handed Cartesian coordinate system strictly aligned with the mechanical characteristics of the physical environment. Simultaneously, when determining the tree height value required for the variance parameter, it does not rely on the yet-to-be-generated digital terrain model, but instead executes a local neighborhood search algorithm: using the current feature point... Centered on the initial sparse point cloud, all points within a horizontal distance of 5.0 meters are searched. The minimum Z-axis coordinate is extracted as the virtual ground elevation. The Z-axis coordinate of the current feature point is subtracted from this virtual ground elevation, and the difference is used as the relative tree height estimate for that feature point. Within a local vertical reference frame constructed with the Z-axis as the reference, an anisotropic variance assignment procedure is executed. The longitudinal variance component along the Z-axis is set to 0.02 m² to respond to the structural stability of vegetation in the direction of gravity, while the lateral variance component perpendicular to the Z-axis plane is set to 1.5 m² to accommodate the horizontal sway error caused by wind load. Then, using the coordinate rotation matrix R formed by the direction cosine matrix of the photogrammetric global coordinate system relative to the local vertical reference frame, the tensor propagation law formula is applied. For anisotropic variance diagonal matrix containing longitudinal and transverse variance components A spatial datum transformation is performed to calculate the anisotropic covariance matrix in off-diagonal form in the photogrammetric global coordinate system. The anisotropic weight matrix P is obtained by inverting the weight matrix P. The anisotropic weight matrix P is substituted into the error equation of the bundle adjustment model as the stochastic model parameter of the observation. In the least squares iterative solution process, the weight matrix constrained adjustment engine releases the horizontal dimension degree of freedom while retaining the tight constraint of the elevation dimension of the vegetation feature points. When the final output exterior orientation elements and ground point three-dimensional coordinates are verified by the checkpoint, it shows that the elevation error of the forest area converges to within 0.15m, and there is no geometric distortion or artifact phenomenon caused by the horizontal swaying of vegetation at the edge of the adjacent rigid building.

[0033] Example 2: To verify the accuracy retention capability and engineering practicality of the anisotropic weighted bundle adjustment model proposed in this invention under non-rigid ground cover interference, a semi-physical simulation test platform containing multiple typical land cover types was constructed. The data source of this platform is based on real-world UAV aerial imagery data (image size 6000×4000 pixels, ground resolution GSD of 0.05m) and the ground truth of the digital surface model (DSM) obtained by high-precision LiDAR scanning. The test scenario is set as a mixed terrain area continuously affected by crosswinds, where the wind speed is set to 8m / s to simulate severe convective weather. This wind load acts on the broad-leaved forest area in the test area, causing the tree canopy feature points to produce random horizontal swaying with an amplitude of 0.5m to 1.2m. At the same time, Gaussian white noise with a mean of 0 and a variance of 0.5 pixels is superimposed on the image grayscale data to simulate electronic thermal noise in the sensor imaging process; An image set containing dynamic non-rigid perturbations was designed with three experimental control groups, each employing a differentiated weighting strategy, to construct a complete chain of evidence. Control group one adopted the traditional isotropic bundle adjustment strategy, assuming all feature points are rigid and setting the weight matrix of all observations to a scalar multiple of the identity matrix, simulating how existing technologies handle non-rigid features. Control group two, as an over-range control group, introduced anisotropic weighting, but set the lateral variance component of vegetation feature points to 0.05 m², verifying the system's response when constraints are too strict and approach the rigid assumption. The experimental group of this invention implemented the technical solution, setting the lateral variance component of vegetation feature points to 1.5 m² and the longitudinal variance component to 0.02 m² based on the hydrodynamic relationship model between wind speed monitoring data and vegetation height. The decision logic for this parameter setting lies in balancing the tolerance for horizontal swaying of vegetation with the preservation of elevation information.

[0034] After the experiment was started, the feature extraction module extracted a total of 152,000 feature points from the image set. In the intermediate data processing stage, for the feature points in the experimental group of this invention that were semantically labeled as vegetation, the system constructed an anisotropic variance diagonal matrix based on the local vertical reference frame calculated from the positioning and orientation data. Key intermediate process data shows that the global covariance matrix after coordinate rotation matrix R transformation... The data exhibits a non-diagonalization characteristic, with the values ​​of its non-diagonal elements increasing dramatically from 0 to the range of 0.4 to 0.8, while the Z-axis component on the diagonal remains low. This data phenomenon reveals the core mechanism of the invention: through the tensor propagation law, the horizontal uncertainty defined in the local coordinate system is accurately mapped and applied to the X and Y components of the global coordinate system, thereby mathematically achieving the directional release of non-rigid errors. The final adjustment results show differences in accuracy. The root mean square error (RMSE-Z) of the checkpoint elevation output of control group one is 0.85m, and obvious stringy geometric artifacts were observed at the eaves of buildings near the forest area, confirming that the horizontal swaying error of the trees contaminated the rigid structure connected to them through the isotropic weight matrix. The RMSE-Z of control group two is 0.72m, and the number of convergence iterations reaches 15. The results showed that the solution oscillations and convergence difficulties were caused by excessively tight constraints. However, the RMSE-Z of the experimental group of this invention converged stably to 0.12m, which is more than 85% higher than the accuracy of the control group. The number of iterations was only 6, the building edge contours were clear, and the orthogonality was well maintained. The gradient sensitivity analysis of the lateral variance parameter showed that when the lateral variance setting value was lower than 0.3m², the elevation accuracy deteriorated sharply and showed obvious nonlinear cutoff characteristics. When the setting value was in the range of 1.0m² to 2.0m², the RMSE-Z remained in the preferred low-level platform area below 0.15m. However, when it exceeded 5.0m², the condition number of the solution system deteriorated, resulting in a decrease in the overall network strength. This confirmed that the value of the lateral variance component is not arbitrarily chosen, but rather there exists an optimal working window determined by the physical amplitude of non-rigid deformation.

[0035] Example 3: In the process of establishing a local vertical reference frame and configuring anisotropic parameters, the following standardized calibration procedure is used to determine the variance component values, thereby eliminating the uncertainty in parameter setting. This procedure defines the physical attribute inputs of non-rigid feature points: for vegetation categories, the system extracts its average tree height H and canopy sparseness S as state variables; for water categories, it extracts its average wave height W of surface waves as a quantitative index. Based on this, the lateral variance component... The calculation is based on the wind load response model in fluid dynamics. The lateral variance component of vegetation feature points is defined as a function of wind speed v, and its calculation formula is as follows: Where k is the aeroelastic coefficient. In the implementation process in the broad-leaved forest coverage area, the aeroelastic coefficient k is set to 0.01, and the real-time wind speed data of the meteorological station in the test area is connected to dynamically calculate the lateral variance setting value of each feature point to achieve adaptive matching to the environmental conditions. For the water feature point, its lateral variance component is directly set to 100m² to remove the geometric constraints of the horizontal dimension.

[0036] Longitudinal variance components The settings follow a gravity stability grading locking logic. For vegetation feature points, based on the rigid support characteristics of tree trunks in the gravity direction, their longitudinal variance is locked at 0.02 m², which corresponds to the inherent noise level of the measurement system. For water feature points, the longitudinal variance is calculated based on the average wave height W, and the calculation formula is as follows: To accommodate elevation uncertainties caused by water surface fluctuations, this calibration procedure was validated under a set of typical conditions including tree heights ranging from 5m to 20m and wind speeds ranging from 2m / s to 12m / s. The results show that the lateral variance component calculated by the model increases parabolically with wind speed, covering the measured amplitude range of tree crown sway. Furthermore, the mean square error of forest elevation in the corresponding adjustment results remains stable within 0.15m across the entire wind speed range, confirming that the procedure can transform environmental disturbances into deterministic physical compensation parameters.

[0037] Example 4: To construct a geomorphic response database supporting anisotropic parameter settings, an offline calibration process based on a controlled wind tunnel environment was established. Standardized sample trees of typical vegetation such as coniferous forests, broad-leaved forests, and shrubs were selected and placed in a hydrodynamic wind tunnel with gradient wind loads ranging from 2 m / s to 20 m / s. A high-frequency laser scanner was used to capture the horizontal displacement vector sequence of tree crown feature points relative to root anchor points in real time at a frequency of 50 Hz. The aeroelastic coefficient k and the nonlinear cutoff threshold of the response function corresponding to different vegetation categories were obtained by least squares fitting regression. Finally, these measured and calibrated physical parameters were encapsulated in the geomorphic attribute lookup table of the airborne processing system as the benchmark data source for subsequent operations to call parameters based on semantic tags.

[0038] During the on-site deployment phase before conducting large-scale aerial photography missions, a baseline calibration procedure for the micro-meteorological environment of the survey area was initiated. This procedure required the deployment of temporary wind speed monitoring nodes or the connection to high-resolution numerical weather prediction data interfaces at typical terrain locations within the survey area. A time-domain cross-correlation analysis was performed using pre-observation data over a period of time and low-frequency attitude disturbances recorded by the airborne inertial measurement unit. This calculated the time delay compensation value for wind field data transmission and the wind speed vertical profile correction factor, which encompasses the terrain roughness of the survey area. The corrected equivalent canopy wind speed was then used as the basis for this calculation. This variable is injected into the aforementioned horizontal variance calculation model, thereby eliminating the systematic bias in parameter estimation caused by spatial sampling mismatch of meteorological data.

[0039] Example 5: To address the potential impact of anisotropic weight matrix computation on system load in large-scale clustered deployment scenarios, this example constructs a computational complexity-aware adaptive resource scheduling procedure, defines a dynamic load assessment model for computational tasks, and monitors in real time the proportion of non-rigid feature points in the image sequence to be processed. ,when When the load exceeds a preset threshold of 30%, the system determines that the current area is a high-load dynamic scene and automatically triggers the high-performance computing mode. In the high-performance computing mode, the system executes a gridded partitioning strategy based on the spatial distribution density of feature points, dividing the survey area into partitions with sides of length [missing information]. The independent subgrids are calculated, and the anisotropy index of feature points within each subgrid is calculated. Its calculation formula is If the subgrid If the value is higher than 1.5 times the average value of the entire survey area, the subgrid is marked as a critical computing unit. The system prioritizes the allocation of GPU computing cores to perform parallel acceleration of the weight matrix inversion process and uses double-precision floating-point format to ensure numerical stability. Conversely, for regular subgrids, the processing is allocated to the CPU thread pool and single-precision floating-point format is used to reduce memory usage.

[0040] In addition, to address potential parameter drift or environmental baseline changes during long-term operation, the system incorporates an online self-checking and dynamic calibration mechanism, extracting a set of control point residual data and calculating its root mean square error. Once discovered If the wind speed response model parameters exhibit a monotonically increasing trend and exceed the safety threshold of 0.2m three times consecutively, then the current wind speed response model parameters are determined. The system has deviated from the optimal state and triggers a parameter reassessment process. It uses the POS data and image matching residuals within the most recent time window to correct the current value of the aeroelastic coefficient k in reverse using the least squares criterion and update it in subsequent adjustment calculations, thereby ensuring the system's continuous adaptability to environmental changes during long-term operations.

[0041] The embodiments of this application have been described above with reference to the accompanying drawings. Unless otherwise specified, the embodiments and features in the embodiments of this application can be combined with each other. This application is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of this application without departing from the spirit of this application and the scope of protection of this invention, and all of these forms are within the protection scope of this application.

Claims

1. A method for analyzing land spatial images based on photogrammetry, characterized in that, Includes the following steps: The system acquires aerial photographic image sequences of the area to be tested, as well as spatiotemporally synchronized positioning and attitude data with the aerial photographic image sequences. It then performs feature extraction and image matching on the aerial photographic image sequences to generate an initial sparse point cloud containing three-dimensional spatial coordinates. Perform semantic segmentation of ground cover categories on aerial photographic image sequences, and divide the feature points in the initial sparse point cloud into a set of rigid feature points and a set of non-rigid feature points based on the segmentation results. For each feature point in the set of non-rigid feature points, establish a local perpendicular reference frame with the direction of the gravity vector as the Z-axis; Within a local vertical reference frame, a longitudinal variance component along the Z-axis and a transverse variance component along a plane perpendicular to the Z-axis are configured for the feature points, wherein the value of the transverse variance component is set to be greater than the value of the longitudinal variance component, thereby constructing an anisotropic variance diagonal matrix. The anisotropic variance diagonal matrix is ​​transformed from the local perpendicular reference system to the photogrammetric global coordinate system using a coordinate rotation matrix. The anisotropic covariance matrix in off-diagonal form is calculated, and the anisotropic weight matrix is ​​obtained by inverting the anisotropic covariance matrix. The anisotropic weight matrix is ​​substituted into the error equation of the bundle adjustment model as the stochastic model parameters of the observations. The least squares iterative adjustment is performed to output the adjusted exterior orientation elements and the three-dimensional coordinates of the ground points.

2. The method for analyzing land spatial images based on photogrammetry according to claim 1, characterized in that, The steps for establishing a local perpendicular reference frame with the direction of the gravity vector as the Z-axis include: extracting the gravity acceleration vector data from the positioning and attitude data, or calculating the fitting plane normal vector of the initial sparse point cloud in the neighborhood of the current feature point, and defining the obtained vector direction as the gravity vector direction; using the gravity vector direction as the Z-axis of the local perpendicular reference frame, and constructing two orthogonal vectors perpendicular to the Z-axis as the X-axis and Y-axis of the local perpendicular reference frame, respectively; calculating the direction cosine matrix of the photogrammetric global coordinate system relative to the local perpendicular reference frame, and determining the direction cosine matrix as the coordinate rotation matrix.

3. The method for analyzing land spatial images based on photogrammetry according to claim 1, characterized in that, The steps of configuring the longitudinal variance component along the Z-axis and the lateral variance component along the plane perpendicular to the Z-axis for feature points specifically include: identifying the specific land cover category of the non-rigid feature points; if the land cover category is vegetation, then setting the longitudinal variance component to a first variance threshold and setting the lateral variance component to a second variance threshold, wherein the second variance threshold is at least 100 times the first variance threshold, so as to reduce the constraint weight on the horizontal displacement of vegetation feature points in the adjustment solution; if the land cover category is water, then setting the longitudinal variance component to the first variance threshold and maintaining the elevation constraint on water surface feature points.

4. The method for analyzing land spatial images based on photogrammetry according to claim 1, characterized in that, The calculation process of the anisotropic covariance matrix satisfies the following mathematical relationship: ,in, R is the anisotropic covariance matrix, and R is the coordinate rotation matrix. It is the transpose of the coordinate rotation matrix. Let P be the anisotropic variance diagonal matrix; the formula for calculating the anisotropic weight matrix P is: .

5. The method for analyzing land spatial images based on photogrammetry according to claim 1, characterized in that, The least squares iterative adjustment solution also includes a weighting step based on the regional geometric rigidity index: calculating the regional geometric rigidity index of the initial sparse point cloud in the local region, which is used to characterize the ratio of the distribution density of rigid feature points to non-rigid feature points in the current adjustment unit; determining the dynamic balance coefficient based on the regional geometric rigidity index; adjusting the weight ratio of the image observation equation and the auxiliary navigation observation equation in the bundle adjustment model using the dynamic balance coefficient; and increasing the weight ratio of positioning and attitude determination data in the adjustment solution when the regional geometric rigidity index is lower than the preset stability threshold by using the dynamic balance coefficient.

6. The method for analyzing land spatial images based on photogrammetry according to claim 5, characterized in that, The calculation logic for the regional geometric stiffness index is as follows: a topological neighborhood is constructed with the currently solved image frame as the center, and the number of rigid feature points within this neighborhood is counted. The number of non-rigid feature points ; Calculate the proportion of rigid points The calculation formula is: ; The proportion of rigid points As a regional geometric stiffness index, among which The value of is positively correlated with the regional geometric stiffness index.

7. The method for analyzing land spatial images based on photogrammetry according to claim 1, characterized in that, The steps for performing least squares iterative adjustment include: calculating the reprojection residual vector of each feature point in a single iteration; calculating the equivalent weight factor using a robust estimation function based on the anisotropic weight matrix and the reprojection residual vector; iteratively updating the anisotropic weight matrix using the equivalent weight factor to reduce the influence of feature points whose reprojection residual vector exceeds three times the mean square error limit on the adjustment results, until the error equation converges.

8. The method for analyzing land spatial images based on photogrammetry according to claim 1, characterized in that, The steps for performing semantic segmentation of land cover categories on aerial photographic image sequences specifically include: inputting the aerial photographic image sequence into a pre-trained semantic segmentation neural network to generate a pixel-level semantic mask; performing spatial projection matching between the pixel-level semantic mask and the initial sparse point cloud to assign a unique semantic label to each 3D point in the initial sparse point cloud; and performing point cloud filtering based on the semantic labels to remove feature points labeled as dynamically moving objects, retain feature points labeled as buildings and roads as rigid feature points, and retain feature points labeled as woodlands and water surfaces as non-rigid feature points.

9. A method for analyzing land spatial images based on photogrammetry according to claim 1, characterized in that, The method also includes: using the exterior orientation elements optimized by least squares iterative adjustment and the three-dimensional coordinates of ground points to perform dense matching processing on the aerial photographic image sequence to generate a high-density point cloud; constructing a digital land surface model based on the high-density point cloud; and using the aerial photographic image sequence to perform texture mapping on the digital land surface model to generate a real-scene three-dimensional model.

Images