A method and system for real-time interpretation and ground feature labeling of unmanned aerial vehicle photogrammetry images

Abstract

The application discloses a kind of unmanned aerial vehicle aerial survey image real-time interpretation and ground feature marking method and system, it is related to real-time interpretation and ground feature marking field, including acquisition positioning module, visual interference module, space-time evolution module and determination processing module, by real-time acquisition multidimensional positioning solution internal data flow, in combination with the clustering and correlation analysis of characteristic curve, determine pose solution jolt event and its time window of occurrence, by analyzing the dense motion vector field between images and feature loss area, identify the causal interference domain that causes jolt, it is projected to physical space in reverse, determine risk analysis area, the multiple risk factors in this area are comprehensively calculated, generate drift risk space-time evolution prediction graph;Space-time alignment and unified mapping are carried out, and the comprehensive marking risk index of each grid is decoded, and the hierarchical response mechanism is started, to obtain the final effectiveness determination and processing of ground feature marking, fundamentally guarantee the spatial credibility and time sequence availability of marking result.

Description

Technical Field

[0001] This invention relates to the field of real-time interpretation and feature annotation, specifically to a method and system for real-time interpretation and feature annotation of UAV aerial survey images. Background Technology

[0002] In existing technologies, drones, due to their high maneuverability and flexibility, are widely used for long-term hovering or low-speed cruising missions in specific areas (such as construction sites, agricultural plots, and disaster sites) to perform real-time video surveillance, ground feature identification, and geographic information system (GIS) annotation. To achieve accurate positioning on consumer or industrial drones without relying on expensive real-time dynamic differential modules, Visual Inertial Odometry (VIO) technology has become the mainstream solution. However, the reliability of the VIO algorithm highly depends on its core assumptions: the static nature and rich texture of the observed scene. In practical applications, drones often face visually challenging environments, such as flying over textureless water surfaces or solid-color walls, encountering large areas of rapidly moving cloud shadows or other dramatic changes in lighting, or the presence of large, dynamic objects in the scene. These situations can violate the assumptions of the VIO algorithm, leading to large-scale loss or mismatch of visual feature points, which in turn causes turbulence or failure of the positioning solution system, ultimately manifesting as sudden jumps or continuous drift in the drone's pose estimation. This drift can accumulate to several meters or more within just a few minutes, causing the geographic feature labels that the operator sees on the screen to have a serious positional deviation on the corresponding geographic map, i.e., the label drift problem. This greatly reduces the spatial reliability and temporal availability of real-time interpretation and labeling tasks.

[0003] To address the aforementioned annotation drift problem, existing technologies have proposed preliminary solutions. One common approach is to simply monitor the number of feature points output by the VIO system; when the number falls below a certain static threshold, the localization is considered unreliable. This single-indicator judgment is too coarse, failing to distinguish between normal feature point fluctuations and systemic failure symptoms, easily leading to numerous false alarms and missed alarms. Another approach attempts to verify the spatial consistency between current features and historically labeled features, but this is a passive and lagging check; when the overall location has already drifted, this internal verification loses its benchmark. More importantly, these methods generally lack the ability to accurately diagnose localization failure events, trace root causes, and proactively predict risks. They cannot immediately initiate a reliable pose compensation strategy to "fill" the data gaps during periods of unreliable visual signals when a failure occurs, forcing annotation work to be interrupted or generating a large amount of invalid data. Therefore, there is a need for a method and system that can perform in-depth state monitoring of VIO systems, accurately diagnose turbulence events, predict future risks, and initiate graded response and pose compensation closed loops without relying on expensive hardware upgrades and through algorithmic optimization, so as to fundamentally solve the pain points of consumer drones in professional surveying and long-term monitoring applications.

[0004] To address the aforementioned shortcomings, a technical solution is provided. Summary of the Invention

[0005] To address the technical problems mentioned in the background section, this invention is proposed. This invention provides a method and system for real-time interpretation and feature annotation of UAV aerial survey images.

[0006] The objective of this invention can be achieved through the following technical solution: a method for real-time interpretation and feature annotation of UAV aerial survey images, comprising the following steps:

[0007] By collecting and analyzing the internal data flow of multidimensional positioning solution in real time, the trajectory drift point is identified using the dynamic baseline model. Combined with the clustering and correlation analysis of characteristic curves, the pose solution turbulence event and its occurrence time window are determined.

[0008] By utilizing the determined turbulence events and their time windows, and by analyzing the dense motion vector field and feature loss areas between images, the causal interference domain of turbulence is identified, and a geospatial distribution map of visual interference sources is generated.

[0009] Based on the geospatial distribution map, the risk analysis area is determined. Taking into account the risk factors of lack of basic texture, changes in lighting, and the occupation of dynamic targets, a spatiotemporal evolution prediction map of drift risk is generated.

[0010] The geospatial distribution map of visual interference sources is spatiotemporally aligned and uniformly mapped with the spatiotemporal evolution prediction map of drift risk. The comprehensive annotation risk index of each grid is decoded to obtain the final validity judgment and processing of the ground feature annotation.

[0011] Furthermore, the steps for determining the pose calculation turbulence event and its occurrence time window are as follows:

[0012] Based on the internal data flow of multidimensional localization calculation, the first visual feature point distribution entropy curve, the second filter state covariance curve, and the third scale factor curve are constructed.

[0013] Based on the analysis of single-source nominal trajectory drift points, a set of discrete drift time points is obtained, and a drift analysis window is defined by clustering.

[0014] Furthermore, the step of determining the pose calculation turbulence event and its occurrence time window also includes:

[0015] Within the drift analysis window, identify the first causal instability characteristic relationship and the second symptom resonance characteristic relationship to determine the pose solution turbulence event. Record the start time and duration of the pose solution turbulence event. The start time is the moment when the distribution entropy curve of the first visual feature point appears to have a significant local minimum. The duration is the time from the start time until the value of the second filter state covariance curve falls back from its significant local maximum to below the upper limit threshold of the determinant value of the fusion filter state covariance submatrix.

[0016] Furthermore, the steps for obtaining the single-source nominal trajectory drift point are as follows:

[0017] When the UAV performs aerial surveying missions, it collects and synchronizes the internal data stream of multidimensional positioning calculation in real time, including the visual odometry scale factor, the determinant value of the fusion filter state covariance submatrix, and the distribution entropy of effective visual feature points.

[0018] Based on the internal data flow of multidimensional positioning solution, dynamic baseline models based on sliding windows are established to identify outliers that exceed a preset threshold and obtain single-source nominal trajectory drift points.

[0019] Furthermore, the steps for identifying the first causal instability characteristic relationship and the second symptom resonance characteristic relationship are as follows:

[0020] Significant local minima were detected on the distribution entropy curve of the first visual feature point, and significant local maxima were detected on the covariance curve of the second filtered state, which was determined to identify the first causal unstable feature relationship.

[0021] Centered on the significant local maximum moment detected in the first causal instability characteristic relationship, the local variance of all data points of the third scale factor curve within the time window of a fixed length is calculated, and it is determined to be the second symptom resonance characteristic relationship.

[0022] Furthermore, the steps for obtaining the geospatial distribution map of the visual interference sources are as follows:

[0023] Full-spectrum filtering is performed on images in high frame rate image sequences before and after the event to generate a set of filtered responses. Based on the response aggregation, a texture intensity energy map is obtained. The texture intensity energy map is then subjected to low-threshold segmentation and morphological optimization to obtain the low-valley region of the full-spectrum filtering response.

[0024] When the spatial overlap rate of the union of the non-dominant motion manifold region and the trough region of the full-spectrum filter response is higher than the preset causal association threshold in the feature loss region, the union is identified as the causal attribution interference domain.

[0025] By obtaining the interference type and influence intensity of the causal attribution interference domain, and by back-projecting it from the image space and aggregating it into the geographic space based on the UAV pose, a geographic spatial distribution map of the visual interference source is obtained.

[0026] Furthermore, the steps for obtaining the non-dominant motion manifold region are as follows:

[0027] The problem time window is determined by solving the turbulence event using the determined pose, so as to obtain the high frame rate image sequence before and after the event occurs;

[0028] The dense motion vector field between two consecutive frames in a high frame rate image sequence before and after the event occurs is calculated. The motion mode is decoupled and the saliency is determined from the dense motion vector field to obtain the non-dominant motion manifold region.

[0029] Furthermore, the steps for obtaining the spatiotemporal evolution prediction map of drift risk are as follows:

[0030] A three-dimensional field-of-view cone is constructed based on a preset task file, and its intersection with the three-dimensional digital surface model of the risk quantification analysis area is calculated in three-dimensional space to obtain the UAV field-of-view footprint.

[0031] Furthermore, the step of obtaining the spatiotemporal evolution prediction map of drift risk also includes:

[0032] After normalizing the illumination change risk index and the dynamic target interference risk index, they are weighted and fused with the basic texture scarcity index to calculate the total risk potential of each spatial location at different times.

[0033] Based on the preset task file, the total risk potential of all locations falling within the UAV's field of view footprint is summed at each moment to obtain the instantaneous risk score at that moment. These scores are then arranged in chronological order to obtain a prediction map of the spatiotemporal evolution of drift risk.

[0034] Furthermore, the steps for obtaining the dynamic target interference risk index are as follows:

[0035] Risk quantification analysis areas were determined based on the geospatial distribution map of visual interference sources;

[0036] The risk quantification analysis area is subjected to grid analysis, the corner density of each grid cell is calculated, and the basic texture scarcity index is obtained by normalizing and reversing the corner density.

[0037] Furthermore, the step of obtaining the dynamic target interference risk index also includes:

[0038] A binary shadow map is generated based on the time-series solar position and a three-dimensional model. The light change risk index is obtained by calculating the time gradient of the binary shadow map.

[0039] A probabilistic dynamic target occupancy grid model is established based on historical data, and the probability and expected velocity of each grid being occupied by a target are predicted to calculate the dynamic target interference risk index.

[0040] Furthermore, the final validity determination and processing steps for the aforementioned feature annotations are as follows:

[0041] By spatiotemporally aligning and uniformly mapping the geospatial distribution map of visual interference sources and the spatiotemporal evolution prediction map of drift risk, a multi-layer defect information map is obtained.

[0042] Based on the multi-layer defect information map, a defect coupling enhancement representation vector is constructed. The defect coupling enhancement representation vector is then used to obtain a comprehensive labeled risk index through risk decoding.

[0043] Based on the comprehensive annotation risk index, a risk classification response and pose correction closed loop is used to make the final validity judgment and processing of ground feature annotations.

[0044] A real-time interpretation and feature annotation system for UAV aerial survey images includes:

[0045] The data acquisition and positioning module collects and analyzes the internal data flow of multidimensional positioning solution in real time, identifies trajectory drift points using a dynamic baseline model, and determines pose solution turbulence events and their time windows by combining clustering and correlation analysis of characteristic curves.

[0046] The visual interference module utilizes the determined turbulence events and their time windows to identify the causal interference domain of turbulence by analyzing the dense motion vector field and feature loss areas between images, and generates a geospatial distribution map of visual interference sources.

[0047] The spatiotemporal evolution module determines the risk analysis area based on the geospatial distribution map, and generates a spatiotemporal evolution prediction map of drift risk by taking into account risk factors such as lack of basic texture, changes in lighting, and the occupation of dynamic targets.

[0048] The judgment and processing module performs spatiotemporal alignment and unified mapping between the geospatial distribution map of visual interference sources and the spatiotemporal evolution prediction map of drift risk, and decodes the comprehensive annotation risk index of each grid to obtain the final validity judgment and processing of the ground feature annotation.

[0049] Compared with the prior art, the beneficial effects of the present invention are:

[0050] This invention acquires and analyzes the internal data stream of multidimensional positioning calculation in real time, identifies trajectory drift points using a dynamic baseline model, and determines pose calculation turbulence events and their time windows by combining clustering and correlation analysis of feature curves. Using the identified turbulence events and their time windows, it backtracks and extracts corresponding high-frame-rate image sequences. By analyzing the dense motion vector field and feature loss regions between images, it identifies the causal interference domain causing the turbulence and projects it back into physical space to generate a geospatial distribution map of visual interference sources. This achieves a paradigm shift from "passive response" to "active prediction and accurate attribution," enhancing the depth and breadth of risk perception. By constructing a multidimensional dynamic baseline model based on visual odometry scale factor, fusion filter covariance, and feature point distribution entropy, this invention can capture real-time "turbulence" events in the positioning system with high accuracy and low false alarm rate. More importantly, this invention innovatively traces back image sequences and, by analyzing dense motion vector fields and full-spectrum filtering responses, can accurately attribute abstract algorithm failure events to specific interference sources in the physical world (such as dynamic objects and weakly textured areas), generating a geospatial distribution map of visual interference sources that can be reused long-term. Combined with quantitative predictions of multiple risks such as future mission route lighting and dynamic targets, this invention elevates the understanding of location risks from post-event discovery to pre-event prediction, providing a solid data foundation and decision-making basis for achieving truly reliable real-time annotation.

[0051] This invention determines risk analysis areas based on the geospatial distribution map of visual interference sources. It comprehensively calculates multiple risk factors within these areas, including lack of basic texture, lighting variations, and dynamic target occupancy. Combined with the UAV's field-of-view footprint along a pre-set mission flight path, it generates a spatiotemporal evolution prediction map of drift risk. The geospatial distribution map of visual interference sources and the drift risk prediction map are spatiotemporally aligned and uniformly mapped to decode the comprehensive annotation risk index for each grid. A tiered response mechanism is then activated to determine the final validity and processing of the ground feature annotations. This constructs an intelligent closed-loop system from risk perception to pose compensation, fundamentally ensuring the spatial reliability and temporal availability of the annotation results. Faced with positioning drift, existing technologies typically only passively mark data as unreliable or interrupt the task, failing to provide effective solutions when problems occur. This invention establishes a closed-loop mechanism for tiered response and pose correction. When the system determines that the risk level exceeds a critical value based on historical and predicted data, it does not simply abandon the process but immediately triggers a "pose smoothing extrapolation" strategy. This strategy seamlessly switches to a high-frequency inertial measurement unit (IMU) for precise integration during a brief window of unreliable visual signals (typically 1-2 seconds) to generate a smooth pose trajectory that conforms to the laws of physical motion to fill data gaps. This proactive, real-time compensation mechanism ensures the continuity and accuracy of the UAV's pose, directly addressing the core pain point of "annotation drift" caused by abrupt pose changes or drift. Therefore, without increasing the cost of expensive hardware (such as RTK), this invention significantly improves the operational quality and data reliability of consumer or industrial UAVs in professional surveying and long-term monitoring applications through deep innovation at the algorithm level, demonstrating practical engineering value and technological originality. Attached Figure Description

[0052] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. The following drawings are not drawn to scale according to the actual size, but are intended to illustrate the main idea of ​​the present invention.

[0053] Figure 1 This is a flowchart of the method of the present invention;

[0054] Figure 2 This is a system block diagram of the present invention;

[0055] Figure 3 This is a flowchart of step S14 of the present invention. Detailed Implementation

[0056] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are also within the scope of protection of the present invention.

[0057] Example 1

[0058] like Figure 1 As shown, a method for real-time interpretation and feature annotation of UAV aerial survey images includes:

[0059] Step S10: By collecting and analyzing the internal data flow of multidimensional positioning solution in real time, the trajectory drift point is identified using the dynamic baseline model, and the pose solution turbulence event and its occurrence time window are determined by combining the clustering and correlation analysis of the characteristic curves.

[0060] Further, step S10 includes:

[0061] Step S11: When the UAV is performing aerial surveying tasks, collect and synchronize the internal data stream of multi-dimensional positioning solution in real time, including visual odometry scale factor, determinant value of fusion filter state covariance submatrix and effective visual feature point distribution entropy.

[0062] The visual odometry scale factor is defined as a key scalar parameter in a monocular visual inertial odometry (VIO) system used to convert the relative scale motion estimated by the vision module into the scale of the real physical world. When a UAV performs aerial surveying missions, it simultaneously collects two types of core data to support the VIO scale factor solution. One type is a continuous image sequence acquired by the onboard industrial camera. This data contains scene texture and ground feature information of the survey area. By extracting feature points and performing matching and tracking, it provides visual motion constraints for the scale factor. The other type is three-dimensional acceleration and three-dimensional angular velocity data acquired by the built-in inertial measurement unit (IMU). This data inherently possesses physical unit attributes, and through integration, the absolute motion increment of the UAV can be obtained, providing an inertial motion reference for the scale factor. The system fuses these two types of data in a loosely or tightly coupled manner to obtain the VIO scale factor. The VIO scale factor solution is a mature existing technology and a common method in current UAV VIO positioning systems; related details will not be elaborated here.

[0063] The fusion-filtered state covariance matrix used by UAVs in aerial surveying missions is a quantitative description of the uncertainty of their output position estimates, provided by algorithms such as the extended Kalman filter. Its function is to directly reflect the reliability of the positioning results. In aerial surveying flights that rely on visual inertial odometry (VIO) for positioning over extended periods, due to the inherent cumulative error of VIO, the values ​​of the position components in this covariance matrix will continuously increase with time and flight distance. This directly indicates that the drift of the positioning results is intensifying, and its reliability is decreasing accordingly. The data acquisition process follows industry standards. For example, within the ROS (Robot Operating System) framework, the node performing sensor fusion publishes messages of type nav_msgs / Odometry to a standard topic. External modules subscribe to these messages and parse the required data from their pose.covariance member (a 36-element array representing a 6×6 covariance matrix). To obtain a single quantitative indicator that comprehensively assesses overall positioning uncertainty, a 3×3 position covariance submatrix specifically characterizing three-dimensional position uncertainty needs to be extracted from the 6×6 covariance matrix, and its determinant value needs to be calculated. This determinant value physically corresponds to the volume of the uncertainty ellipsoid, integrating the variances of the X, Y, and Z axes and their covariances, thus serving as a comprehensive scalar value for real-time positioning quality assessment. The step of obtaining the determinant value of the fusion filter state covariance submatrix is ​​a mature existing technology and a common method in current UAV VIO positioning systems; related details will not be elaborated here.

[0064] The current image frame from the real-time video stream being processed by the Visual Odometry (VIO) system front-end is acquired. This image frame is processed to extract the pixel coordinates of all successfully tracked valid visual feature points. Based on the distribution of these coordinate points on the image plane, a quantized value, namely the valid visual feature point distribution entropy, is obtained through spatial discretization and information entropy calculation. Its function is to evaluate the quality of visual information used for localization in the current frame in real time. This method is a mature existing technology in the field of visual odometry (VIO) and will not be elaborated on here. The visual odometry scale factor, the fusion filter state covariance matrix, and the valid visual feature point distribution entropy are synchronously acquired and accompanied by precise timestamps, forming the internal data stream for multidimensional localization calculation.

[0065] Based on the internal data flow of multidimensional localization calculation, dynamic baseline models based on sliding windows are established to identify outliers exceeding a preset threshold, thus obtaining single-source nominal trajectory drift points.

[0066] Specifically, in step S12: Based on the internal data flow of multi-dimensional positioning solution, establish dynamic baseline models based on sliding windows. By calculating the Z-score of the current value and the dynamic baseline model, identify outliers that exceed the preset threshold and obtain the single-source nominal trajectory drift points.

[0067] Step S12 processes the three independent time-series data streams of the multidimensional localization solution internal data stream to identify potential system instability symptoms caused by a single data source. The core of this step is to establish a dynamic baseline model for each time-series data stream that can adapt to the current system operating conditions. The reason for using a dynamic baseline instead of a static threshold is that the normal fluctuation range of various internal indicators of the VIO system changes continuously under different motion states (e.g., stationary, uniform speed, acceleration / deceleration) and different environments (e.g., rich texture, sparse texture). A dynamic baseline can more accurately reflect the normal level under the current conditions. The process of establishing this sliding window-based dynamic baseline model is as follows: the system maintains a fixed-length FIFO data queue of length N as a sliding window, for example, N=100, to store the values ​​of this indicator for the most recent N moments. At each new moment t, when a new data point X... t When data is collected, the system calculates the arithmetic mean μ of all data points within the current window. t and standard deviation σ t Arithmetic mean μ t and standard deviation σ t Together, they form the dynamic baseline model at time t. As time progresses, the window slides forward, discarding old data and incorporating new data, μ t and σ t This is also updated in real time, allowing the baseline model to closely follow the short-term trends and volatility changes of the data sequence, dynamically defining the normal range. After establishing the dynamic baseline model, the system calculates the current value X. t The deviation Z is quantified by the Z-score (standardized score) of this dynamic baseline. t The calculation formula is Z. t =(X t -μ t ) / σ t The Z-score is calculated by comparing the distance of a current data point from its recent mean to its recent standard deviation. The absolute value of the calculated Z-score is then compared to a preset threshold to identify statistically significant outliers. Here, a threshold of 3 is chosen. The preset threshold is based on empirical rules in statistics, such as the 68-95-99.7 rule. In a dataset that approximates a normal distribution, approximately 99.7% of the data points will fall within 3 standard deviations of the mean (i.e., ...). Within the range of ), the absolute value of the Z-score of a data point exceeds 3. This means it's an extremely low-probability event with a probability of less than 0.3% occurring in the current dynamic environment. (Option 3) As a threshold, it achieves a good balance between detection sensitivity and false alarm rate: it is strict enough to effectively filter out most normal random noise and slight fluctuations, avoiding frequent false alarms caused by an excessively low threshold; yet it is sensitive enough to promptly capture extreme deviations that truly indicate significant changes in the system state, avoiding missed detections caused by an excessively high threshold. When the system identifies that the absolute value of the Z-score of a time series at a certain moment exceeds 3... The threshold is then used to mark the data point at that moment as a single-source nominal trajectory drift point. By performing the above process on the three dimensions of scale factor, location variance, and distribution entropy respectively, the system can independently identify single-source nominal trajectory drift points in each dimension, providing a basis for subsequent comprehensive system status assessment and fault diagnosis.

[0068] Step S13: Construct the first visual feature point distribution entropy curve, the second filtered state covariance curve, and the third scale factor curve based on the internal data flow of the multidimensional localization solution;

[0069] Based on the analysis of single-source nominal trajectory drift points, a set of discrete drift time points is obtained, and a drift analysis window is defined by clustering.

[0070] During UAV aerial surveying missions, the internal data stream output in real time by the multi-dimensional positioning and calculation module is monitored. A time-sliding window method is employed to continuously sample and record multiple key performance indicators within this data stream, constructing trend curves reflecting the dynamic changes in system performance from discrete instantaneous data points. Specifically, the real-time acquired effective visual feature point distribution entropy, the determinant of the fusion filter state covariance submatrix, and the visual odometry scale factor are used as data points in the time series. Based on the continuous sampling values ​​within this time-sliding window, the first visual feature point distribution entropy curve, the second filter state covariance curve, and the third scale factor curve are constructed, respectively.

[0071] The discrete set of drift time points is obtained based on the drift point analysis of the single-source nominal trajectory. Clustering is used to define drift analysis windows with analytical value. This process aims to merge multiple discrete time points that are temporally adjacent and belong to the same physical event into a single independent "bump event." A time threshold is set. (For example (seconds), this value defines the standard for judging whether time points are close together, for All time points in the list Sort the events in chronological order, create an empty list of event clusters, and iterate through all subsequent event clusters starting from the first event. And calculate its relationship with the previous time point. Time difference ,like This indicates the time point. and Belonging to the same event, If added to the current event cluster, This indicates the start of a new event, the end of the current event cluster construction, and the beginning of a new event cluster. Create a new event cluster as the starting point. Repeat this process until all time points have been assigned to the corresponding event clusters. The range of time points contained in each formed event cluster is used to obtain the drift analysis window.

[0072] like Figure 3 As shown, step S14: Identify the first causal instability feature relationship and the second symptom resonance feature relationship within the drift analysis window to determine the pose solution turbulence event, and record the start time and duration of the pose solution turbulence event. The start time is the moment when the distribution entropy curve of the first visual feature point appears to have a significant local minimum value, and its duration is the time from the start time point until the value of the second filter state covariance curve falls back from its significant local maximum value to below the upper limit threshold of the determinant value of the fusion filter state covariance submatrix.

[0073] Identifying the primary causal instability characteristic relationship and the secondary symptom resonance characteristic relationship also includes the following steps:

[0074] Step S141: A significant local minimum is detected on the distribution entropy curve of the first visual feature point, and a significant local maximum is detected on the covariance curve of the second filtered state, which is determined to identify the first causal unstable feature relationship.

[0075] Centered on the significant local maximum moment detected in the first causal instability characteristic relationship, the local variance of all data points of the third scale factor curve within the time window of a fixed length is calculated and determined to be the second symptom resonance characteristic relationship.

[0076] Specifically, step S142: Taking the moment of the significant local maximum detected in the first causal instability characteristic relationship as the center, calculate the local variance of all data points of the third scale factor curve within the time window through a fixed-length time window. If the local variance value is significantly higher than its historical average level, it is determined to be the second symptom resonance characteristic relationship.

[0077] To ultimately determine a pose calculation turbulence event and its temporal characteristics, all the following identification processes must be performed within the drift analysis window. Within this window, the first causal instability characteristic relationship is identified. The logic is that deterioration in sensor information quality (decreased entropy) leads to increased system positioning uncertainty (increased covariance). Its morphological characteristics are: a significant local minimum is detected on the first visual feature point distribution entropy curve, and simultaneously, a significant local maximum is detected on the second filtered state covariance curve. A data point is considered a significant local minimum if two conditions are met: 1) the value of this point must be lower than the values ​​of its two immediately preceding and following data points on the time axis, thus forming a local low point; 2) the value of this point itself must be lower than a pre-set lower limit threshold for the visual feature point distribution entropy, representing severe deterioration in information quality. A data point is considered a significant local maximum if it meets two conditions: 1) the value of this point must be higher than the values ​​of its two immediately preceding and following data points on the time axis, forming a local high; 2) the value of this point must exceed a pre-set upper limit threshold for the determinant of the fusion filter state covariance submatrix, representing excessively high system uncertainty. The moments of significant local minima on the entropy curve are recorded. Significant local maxima of the sum and covariance curve Subsequently, the condition for determining causal instability is that these two moments must be closely coupled in time, i.e., satisfying the condition... ,in The presence of this feature pair, representing a minimal time tolerance for system response delay, is considered the first strong evidence of a turbulent event. Identifying the second symptom resonance feature relationship within the drift analysis window is based on the synchronous manifestation of extreme uncertainty within the system and severe jitter (scale factor fluctuation) in the solution results. This is achieved by using the significant local maxima detected in the first causal instability feature relationship. Define a fixed-length time window centered on [the element]. Calculate the third scale factor curve within this time window. The local variance of all data points within a window is a standard statistic measuring the dispersion of the scale factor data within that window. The condition for symptom resonance is that the calculated local variance must be significantly higher than its historical average. This is determined by a quantification rule: comparing the local variance of the current window with a historical average variance obtained during normal system operation; the local variance is considered significantly higher if and only if the current local variance is greater than a certain percentage of the historical average variance. When doubled, among which For example, the adjustable sensitivity coefficient. The occurrence of a pose resolution turbulence event is considered the second strong piece of evidence if the criteria are met. A pose resolution turbulence event can be determined with high confidence if and only if both causal instability and symptom resonance are observed simultaneously within the same drift analysis window. The start time of this event is defined as the moment when the entropy curve of the first visual feature point distribution shows a significant local minimum. Its duration is the time from the starting point until the value of the second filter state covariance curve falls from its significant local maximum to below the upper limit threshold of the determinant value of the fusion filter state covariance submatrix.

[0078] Steps S14, S141, and S142 together constitute the core of accurate judgment of pose calculation turbulence events. Their analysis is based on establishing a rigorous logical chain from cause to effect to symptom, identifying with high confidence the systemic problems truly affecting positioning accuracy. Step S141 is based on the "causal instability characteristic relationship," meaning that the deterioration of visual information quality (cause: significant local minima appearing in the feature point distribution entropy) inevitably leads to a sharp increase in system positioning uncertainty (effect: significant local maxima appearing in the fusion filter state covariance). This causal correlation analysis based on physical meaning avoids misjudgments caused by single indicator anomalies, constituting the first strong piece of evidence for event judgment. Step S142 is based on the "symptom resonance characteristic relationship," meaning that when the system is extremely uncertain (covariance reaches its peak), its external manifestation will inevitably be severe jitter in the calculation results (significantly increased local variance of the scale factor). This is equivalent to verifying the instability of the internal state through external manifestations, constituting the second strong piece of evidence for event judgment. Step S14 combines the previous two steps, requiring that both "causal instability" and "symptom resonance" conditions be met simultaneously within the same drift analysis window to ultimately determine it as a "pose calculation turbulence event" and precisely define its occurrence time window (start time and duration). This combination of steps is indispensable to the entire method because it plays a crucial role in "problem discovery and definition." Without this series of steps, subsequent processes will not be able to obtain a clear and high-confidence "problem time window." This will directly prevent step S20 from starting, meaning it will be impossible to trace back and extract the high frame rate image sequence corresponding to the problem period for in-depth root cause analysis. The entire system will degenerate into a rudimentary system that can only vaguely perceive drift but cannot pinpoint the exact moment the problem occurred, let alone trace the root cause. Subsequent steps such as visual interference source identification (S20), risk prediction (S30), and final risk classification response (S40) will be impossible due to the lack of accurate "event" input, and the effectiveness of the entire method will cease to exist.

[0079] A method for real-time interpretation and feature annotation of UAV aerial survey images also includes:

[0080] By utilizing the determined turbulence events and their time windows, and by analyzing the dense motion vector field and feature loss areas between images, the causal interference domain of turbulence is identified, and a geospatial distribution map of visual interference sources is generated.

[0081] Specifically, in step S20: using the determined turbulence events and their time windows, backtrack and extract the corresponding high frame rate image sequences, and by analyzing the dense motion vector field and feature loss areas between images, identify the causal interference domain that caused the turbulence, and project it back into the physical space to generate a geospatial distribution map of the visual interference source.

[0082] Specifically, step S20 also includes:

[0083] Step S21: Use the determined pose to solve the turbulence event to determine the problem time window, so as to obtain a high frame rate image sequence before and after the event occurs.

[0084] While calculating the effective visual feature point distribution entropy for each frame in real time, the system actually records this high frame rate image along with its precise capture timestamp in an internal data log that can be queried retrospectively. When the system monitors a turbulence event during pose resolution, it captures and records the precise start and end times of the event, forming a clear problem time window. Using this time window as a query index, the system retrieves and extracts all consecutive image frames whose timestamps fall within this window from the data log containing all historical images, thus locating and obtaining a high frame rate image sequence covering the entire process before and after the turbulence event.

[0085] Step S22: Calculate the dense motion vector field between two consecutive frames in the high frame rate image sequence before and after the event, and perform motion mode decoupling and saliency discrimination on the dense motion vector field to obtain the non-dominant motion manifold region.

[0086] Multi-scale, multi-directional full-spectrum filtering is performed on images in high frame rate image sequences before and after the event to generate a set of filter responses. Based on the response aggregation calculation, a texture intensity energy map is obtained. The texture intensity energy map is then subjected to low-threshold segmentation and morphological optimization to obtain the low-valley region of the full-spectrum filter response.

[0087] Two consecutive frames are extracted chronologically from a high frame rate image sequence before and after the event, denoted as . and The reason for choosing consecutive frames is to capture the most instantaneous, pixel-level motion in the scene. Dense optical flow algorithms, such as the classic Gunnar Farnebäck algorithm, are used to analyze these two frames. For each pixel in the array, calculate its position. At the corresponding position in the graph, a two-dimensional motion vector is generated for each pixel. This vector describes the pixel's position in the horizontal direction between frames. and vertical direction The displacement on the image yields a dense motion vector field with the same dimensions as the original image. This is a two-dimensional data set with the exact same size as the original image. Each element in the set corresponds to a two-dimensional motion vector of a pixel in the original image. To extract meaningful motion patterns from this chaotic vector field, an unsupervised clustering algorithm, such as K-Means, is used. This algorithm clusters the motion vectors of each pixel... As a feature, all pixels in the image are grouped, with pixels corresponding to motion vectors of similar direction (angle) and size (magnitude) being grouped into the same cluster. Background pixel movement caused by camera motion forms a dominant motion cluster occupying a large portion of the image area, i.e., the background motion cluster. The average motion vector of this background motion cluster is calculated and denoted as... Traverse all other non-background motion clusters. And calculate their respective average motion vectors. Motion pattern decoupling refers to separating the dominant background and non-dominant foreground motion patterns that are mixed together. A cluster To be classified as abnormal motion, it must simultaneously meet the significance criteria, including two quantitative conditions: 1) its motion pattern differs significantly from the background motion, as determined by calculating the Euclidean distance between the two mean vectors. To measure, that is ,in 1) The preset motion difference threshold; 2) The cluster must have a sufficient size to exclude the influence of random noise or small objects, the criterion being the total number of pixels occupied by the cluster. With the total number of pixels in the image The proportion exceeds a preset area threshold. ,Right now ,For example It can be set to 0.01, representing 1% of the total image area. The regions covered by all pixel clusters that simultaneously satisfy both conditions are merged to form the non-dominant motion manifold region.

[0088] This is performed on single frames within a high frame rate image sequence. For example, it can be processed through a pre-built Gabor filter bank, a powerful texture analysis toolbox containing a large number of filters for different orientations. and different scales Gabor filters (or spatial frequency) For example, the orientation can take values ​​in 30° increments from 0° to 180°, while setting multiple different scales to capture a variety of textures from coarse to fine. Input image With each filter in the filter bank Perform two-dimensional convolution operations ( Each convolution generates a response map. The value of each pixel in the image This represents the location of the original image. The response intensity of a pixel to the texture features it is sensitive to for that specific filter. A region lacking texture is essentially one that lacks significant edges or patterns in any direction and at any scale. In such regions, the response values ​​of all Gabor filters will be very low. To synthesize information from all filters, the response values ​​for each pixel location in the image are calculated. Calculate a total energy value to generate a texture energy map. This energy value By taking this pixel The response is calculated using the maximum absolute value of the response across all response graphs. Or by calculating the sum of squares of the responses, Where max represents the maximum value. This energy map visually reflects the texture richness across the image. Setting a lower energy threshold... All in the energy diagram Pixels with energy values ​​below this threshold satisfy the condition. The pixels were initially identified as having missing textures. After morphological post-processing, the region formed by these pixels was identified as a low-level area of ​​the full-spectrum filter response. Morphological post-processing included filling internal holes with closing operations and removing isolated noise points with opening operations.

[0089] Step S23: When the spatial overlap rate of the union of the non-dominant motion manifold region and the trough region of the full-spectrum filter response is higher than the preset causal association threshold in the feature loss region, the union is identified as the causal attribution interference domain.

[0090] Based on the high frame rate image sequence before and after the event, the frame preceding the turbulence event is identified and denoted as frame [frame number missing]. From visual inertial odometry Extract the set of valid 3D map points that are currently being successfully tracked from the system. These points were selected based on the following criteria: their reprojection error in the current frame is less than 2.0 pixels, and they have been successfully tracked for more than 5 consecutive frames. A fixed camera intrinsic parameter matrix, obtained in advance through the camera calibration process, is used. and based on The system's motion estimation algorithm, such as extended Kalman filtering, outputs data in real time from frames. Up to the current frame Relative pose transformation matrix , will set Each three-dimensional point in Projected onto the current frame On the image plane, a two-dimensional projection point set is obtained. For each point in this projection point set, in the current frame... With the next frame Optical flow tracing is performed between frames, and a forward-reverse consistency check is used to determine if the tracing has failed. Specifically, if a point from frame ... Tracked to frame Then reverse track back the frame If the Euclidean distance between the current location and its original location is greater than 1.0 pixels, then the point is considered a tracking failure. The convex hull formed by all these tracking failure points on the current frame image is calculated to form the feature loss region. Calculate the region where this feature is lost. Union of non-dominant motion manifold regions and the trough region of full-spectrum filter response Spatial overlap rate between , Indicates the non-dominant motion manifold region. The formula for calculating the spatial overlap rate, representing the trough region of the full-spectrum filter response, is as follows: ,in Represents the pixel area of ​​the region. The calculated overlap rate... Compare with a preset causal association threshold, for example .like If the value exceeds this threshold, it indicates a strong spatial correlation between the large-scale loss of feature points and these two types of interference regions, confirming that the loss is caused by... and The area formed by these two entities constitutes the causal attribution interference domain for this turbulence event.

[0091] Step S24: By obtaining the interference type and influence intensity of the causal attribution interference domain, and based on the UAV pose, back-projecting it from the image space and aggregating it into the geographic space, a geographic spatial distribution map of the visual interference source is obtained.

[0092] To identify the type of interference in the causal attribution interference domain, such as the non-dominant motion manifold region. Or the low point of the full-spectrum filter response Image boundaries and spatial overlap The intensity of the quantified impact, for example, if the calculated for The initial value of the impact intensity of this event. That is Based on the six-degree-of-freedom pose recorded by the drone at that time , This represents the three-dimensional position coordinates of the drone in space. The six parameters—roll, pitch, and yaw attitude angles—represent the UAV's rotation around each axis, together forming the six-degree-of-freedom pose data used for spatial positioning and projection calculations. Through backprojection transformation, this data is precisely mapped from a two-dimensional image coordinate system to a three-dimensional GIS geographic coordinate system, yielding geographic information primitives with geographic boundaries, interference types, and quantization intensity values. When interference sources are repeatedly detected at the same geographic location in multiple missions, these primitives from different times are aggregated. For example, if a region is detected three times consecutively, their influence intensities are... , , The final strength after polymerization The maximum value strategy can be used to update it. This is used to characterize the highest risk level of the area. Within the same geographical location, such as the same GIS grid cell, by performing the above projection and aggregation operations on all detected interference sources, a geospatial distribution map of visual interference sources is obtained. This map, based on a geographic map, shows the visual positioning interference risk level and its attribution type at different geographical locations, providing a detailed minefield map for the UAV's visual navigation system. The back projection transformation, which accurately maps the two-dimensional image coordinate system to the three-dimensional GIS geographic coordinate system, is a standard and commonly used technique in computer vision and photogrammetry, and will not be elaborated here.

[0093] Step S24 is the bridge to achieve the crucial transformation from "image space" to "geospace." It relies on precisely linking the abstract interference analysis results at the image level with real-world geographical locations, thereby giving risk information practical application value. The core of this step is to utilize the precise six-DOF pose data recorded by the UAV to map the type of the "causal attribution interference domain" identified in step S23 and the influence intensity quantified by spatial overlap rate from the two-dimensional image coordinate system to the three-dimensional Geographic Information System (GIS) coordinate system through back projection transformation. Furthermore, it includes an important aggregation processing mechanism: updating the intensity (e.g., taking the maximum value) of multiple interferences detected at the same geographical location at different times, thus establishing a dynamically updated and continuously improving risk database. This step is indispensable in the entire method because it elevates the process from "instantaneous problem diagnosis" to "long-term risk documentation." Without step S24, the "causal attribution interference domain" analyzed in steps S20 to S23 would only exist in a single image or a few images from a single flight, becoming a collection of isolated and unusable data. The subsequent step S30 will be impossible because it requires a "geospatial distribution map of visual interference sources" to determine the risk analysis area and make subsequent risk predictions. Without the "minefield map" generated in step S24, risk prediction becomes unfounded. Ultimately, the entire system will be unable to provide any risk avoidance suggestions based on historical experience for future route planning, nor will it be able to build a global risk situation awareness capability.

[0094] A method for real-time interpretation and feature annotation of UAV aerial survey images also includes:

[0095] Based on the geospatial distribution map, the risk analysis area is determined. Taking into account the risk factors of lack of basic texture, changes in lighting, and the occupation of dynamic targets, a spatiotemporal evolution prediction map of drift risk is generated.

[0096] Specifically, step S30: Based on the geospatial distribution map of visual interference sources, determine the risk analysis area, comprehensively calculate multiple risk factors such as the lack of basic texture, changes in illumination, and the occupation of dynamic targets within the risk analysis area, and combine the UAV field of view footprint of the preset mission route to generate a drift risk spatiotemporal evolution prediction map.

[0097] Specifically, step S30 also includes:

[0098] Step 31: Determine the risk quantification analysis area based on the geospatial distribution map of visual interference sources;

[0099] The risk quantification analysis area is subjected to grid analysis, the corner density of each grid cell is calculated, and the basic texture scarcity index is obtained by normalizing and reversing the corner density.

[0100] A binary shadow map is generated based on the time-series solar position and a three-dimensional model. The light change risk index is obtained by calculating the time gradient of the binary shadow map.

[0101] Based on historical data, a probabilistic dynamic target occupancy grid model is established by region and time period, and the probability and velocity expectation of each grid being occupied by a target are predicted in order to calculate the dynamic target interference risk index.

[0102] The geospatial distribution map of visual interference sources has defined the boundaries for the aerial survey mission and determined the risk quantification analysis area. Static low-texture risk, ambient lighting dynamic shadows, and dynamic interference target analysis are conducted based on this risk quantification analysis area, and are limited to this area. The core of quantifying static low-texture risk lies in calculating the basic texture scarcity index for each spatial location within the risk quantification analysis area. A three-dimensional digital surface model (DSM) of the risk quantification analysis area is constructed using UAV photogrammetry or LiDAR scanning technology. The construction of the DSM is a commonly used existing technique and will not be elaborated upon here. The three-dimensional digital surface model of the risk quantification analysis area is then meshed, for example, into 0.5m × 0.5m grid cells. ,in For grid indexing, the FAST corner density within each grid cell is calculated using high-resolution orthophotos or previous survey images. Specifically, the density of corner points falling within each grid cell is statistically analyzed. Number of FAST corner points within the projection area and divide by the area of ​​the unit. Thus, the corner density is obtained. The density value is then transformed into a standardized index representing the degree of scarcity, and the maximum corner density value is found within the risk quantification analysis region. The density of each unit is normalized to obtain the normalized texture richness. Its range is Basic texture scarcity index For example, an area with extremely sparse texture, such as the surface of calm water, Approaching 0, resulting in Approaching 0, its It then approaches 1.0, indicating extremely high static risk.

[0103] The aerial survey mission date, time period, and geographical latitude and longitude of the risk quantification analysis area are obtained. Based on astronomical algorithms, the data for each second within the mission time period is calculated. Sun azimuth and elevation angle A three-dimensional digital surface model of the risk quantification analysis area was obtained using shadow volume technology, specifically based on the solar azimuth angle. and elevation angle For each occluder in the 3D model, at any given time Rays are emitted from the sun's position to all vertices of the obscuring object, resulting in an infinitely extending shadow volume. By intersecting this shadow volume with a three-dimensional digital surface model, the time of intersection of the ground and building surfaces is calculated. The shadow coverage area is used to obtain a binary shadow map. ,in For grid index, the value is This indicates that the cell is in shadow, and its value is... This indicates that the location is under illumination. Since the risks of visual odometry (VIO) exist not only within shadows but also at the edges where light and shadow rapidly change, the risk index for changes in illumination is... Defined as the time gradient of the shadow map, its specific calculation formula is as follows: , where Δt represents the time difference. This model can accurately predict when the shadow edge will sweep across any specific location within the risk quantification analysis area. This allows for the quantification of navigation risks caused by sudden changes in illumination.

[0104] A probabilistic dynamic target occupancy grid model is established. Specifically, the risk quantification analysis area is divided into different functional zones such as roads, sidewalks, construction sites, and parking lots. Based on historical data, such as urban traffic flow data and surveillance video analysis, statistical models are established for the frequency and average speed distribution of dynamic targets in different functional zones at different times, such as weekday morning rush hour, midday, and nighttime. Based on this model, the risk quantification analysis area is gridded on the ground, and the value of each grid is calculated. At any time in the future The probability of being occupied by a dynamic target Risk is positively correlated with the expected value of the probability of target occupancy and its movement speed; a dynamic target interference risk index is defined. Its calculation formula is ,in It refers to the grid At the moment Expected speed value It refers to the grid At the moment The speed.

[0105] Step 32: Construct a three-dimensional field of view cone based on the preset task file, and perform three-dimensional space intersection calculations between it and the three-dimensional digital surface model of the risk quantification analysis area to obtain the UAV field of view footprint.

[0106] The input consists of a pre-defined mission file containing 3D waypoints, flight speed, camera intrinsic and extrinsic parameters, and camera gimbal attitude control laws. B-spline interpolation is performed on the discrete waypoints in the mission file, and high-frequency sampling is used to obtain each fine time step within the future mission cycle. Six-DOF pose of the unmanned aerial vehicle This pose is determined by three-dimensional spatial coordinates. and three-dimensional rotation angle constitute, These represent the UAV's eastward, northward, and skyward positions in the world coordinate system. The skyward position refers to the U-axis in the ENU coordinate system, which is perpendicular to the ground plane and pointing upwards. These represent the roll, pitch, and yaw angles of the UAV around its own coordinate system's longitudinal (forward), transverse (lateral), and vertical axes, respectively. (The image shows the camera positioned in...) The three-dimensional field-of-view cone at any given moment is defined by using camera intrinsic parameters (such as focal length, principal point, and sensor size) to define a square pyramid in the camera coordinate system. The vertex of this pyramid is the camera's optical center, and the base is the camera's imaging plane. Then, the UAV's pose is used... The camera's mounting attitude relative to the drone's body is considered, and this quadrangular pyramid is transformed from the camera coordinate system to the world coordinate system. The construction of this field-of-view cone is a standard technique in computer graphics and photogrammetry, and will not be elaborated upon here. A geometric intersection operation is performed to obtain the accurate field-of-view footprint. This operation calculates the intersection of the 3D field-of-view cone with a pre-constructed 3D digital surface model of the risk quantification analysis area in 3D space. The intersection formed on the 3D digital surface model is an irregular polygonal region, which represents the area within the... The field of view of a drone camera at any given moment is the area covered by the projection of its view onto various terrain features, i.e., the field of view footprint at that moment. .

[0107] Step 33: After normalizing the illumination change risk index and the dynamic target interference risk index, the total risk potential of each spatial location at different times is calculated by weighted fusion with the basic texture scarcity index. Based on the preset task file, the total risk potential of all locations falling within the UAV's field of view footprint is summed at each time to obtain the instantaneous risk score at that time. The scores are arranged in chronological order to obtain the drift risk spatiotemporal evolution prediction map.

[0108] Risk index for changes in light intensity Risk index of interference with dynamic targets Perform maximum and minimum value normalization processing to map its numerical range to... Intervals are defined to ensure consistency in dimensionality among different risk sources. For any ground grid unit in space... At any time Total risk By analyzing the basic texture scarcity index Normalized light change risk index and dynamic target interference risk index The calculation is performed using a weighted fusion method, and its specific form is as follows: The weighting coefficients among them. The sensitivity of the adopted visual inertial odometry (VIO) calculation method is precisely determined through sensitivity calibration. This calibration process involves independently testing the positioning drift of the VIO system caused by a single risk source (such as weak texture, drastic lighting changes, or dynamic object interference) in a controlled experimental environment. The weights of each risk index are then adjusted based on the magnitude of the drift, ultimately ensuring that the sum of the weight coefficients is 1. For example, if a VIO algorithm is most sensitive to texture, second most sensitive to lighting, and least sensitive to dynamic objects, then the weights can be specifically set as follows: Based on the preset task file, for each fine time step on the preset route. Combined with the drone's field of view footprint at that moment By analyzing all ground grid cells falling within this footprint range Total risk Summing is performed to calculate the instantaneous risk score at that moment, i.e. This score quantifies the VIO positioning drift risk that the image captured by the camera at that moment might cause. All time steps... Calculated The values ​​are arranged in chronological order to obtain a spatiotemporal evolution prediction map of drift risk. This prediction map can be intuitively presented as a one-dimensional curve, with the horizontal axis representing time or flight distance and the vertical axis representing specific data. Numerical values ​​can also be rendered as colored paths along the flight route on a 3D map, with the shades of the path color intuitively reflecting the flight path. The level of risk can be determined to provide sufficient and quantifiable risk warning information for task planning and online decision-making.

[0109] Step S33 is the core of the fusion calculation that realizes the transformation from "multi-source heterogeneous risk" to "unified quantitative prediction." It constructs a comprehensive risk assessment model that can weightedly fuse multi-dimensional risk factors such as static, dynamic, and environmental factors, and combine them with the specific flight path of the future mission route to generate a quantifiable risk prediction curve with spatiotemporal dimensions. This step unifies the dimensions of different risk indices (illuminance, dynamic targets) through normalization, and then, based on the sensitivity calibration results of different VIO algorithms for each risk source, sets weight coefficients for weighted summation to obtain the "total risk potential" of each geographic grid at any given time. Most importantly, it uses the "field of view footprint" of the preset mission as a dynamic sampling window, integrating and summing the total risk potential within the area "seen" by the UAV at each moment to obtain an "instantaneous risk score" that changes over time, forming a "drift risk spatiotemporal evolution prediction map." This step is indispensable in the method because it plays the role of a "risk forecaster," transforming the static risk map into a dynamic risk evolution process closely related to the specific mission. Without step S33, the intermediate data calculated in steps S31 and S32—such as the lack of basic texture, lighting variations, dynamic target risks, and field-of-view footprints—would be fragmented and disjointed, failing to form a unified, future-oriented risk metric. This would prevent subsequent steps S40 (especially S41 and S42) from obtaining the crucial input of the "drift risk defect layer." Without this spatiotemporal evolution prediction map, the system cannot predict the risks of future annotation tasks, thus failing to generate a comprehensive annotation risk index. The hierarchical response mechanism would become passive and lagging due to the lack of forward-looking risk assessment.

[0110] A method for real-time interpretation and feature annotation of UAV aerial survey images also includes:

[0111] The geospatial distribution map of visual interference sources is spatiotemporally aligned and uniformly mapped with the spatiotemporal evolution prediction map of drift risk. The comprehensive annotation risk index of each grid is decoded to obtain the final validity judgment and processing of the ground feature annotation.

[0112] Specifically, step S40: The geospatial distribution map of visual interference sources and the spatiotemporal evolution prediction map of drift risk are spatiotemporally aligned and uniformly mapped to decode the comprehensive annotation risk index of each grid and activate the hierarchical response mechanism to obtain the final validity determination and processing of the ground feature annotations.

[0113] Specifically, step S40 also includes the following steps:

[0114] Step S41: Spatiotemporally align and uniformly map the geospatial distribution map of visual interference sources and the spatiotemporal evolution prediction map of drift risk to obtain a multi-layer defect information map;

[0115] The influence intensity of visual interference sources in the geospatial distribution map is used as the quantification value of visual interference defects. The maximum instantaneous risk score of each geographic grid in the drift risk spatiotemporal evolution prediction map within a preset time period is extracted and used as the quantification value of risk defects. These two types of defect data are spatiotemporally aligned and uniformly mapped to a high-resolution geographic grid coordinate system, resulting in a multi-layered defect information map composed of a visual interference defect layer and a drift risk defect layer. The preset time period is determined based on the characteristic duration of the core task, such as the average execution time of a single aerial survey mission. By setting the preset time period to a specific multiple of the mission duration (e.g., twice), the time range of risk prediction can be ensured to match the needs of actual operations. The purpose of this setting is to assess and quantify the maximum risk that may be encountered in one or more future mission execution cycles, thereby providing targeted, future-oriented risk predictions for decision-making.

[0116] Based on the multi-layer defect information map, a defect coupling enhancement representation vector is constructed. The defect coupling enhancement representation vector is then used to obtain a comprehensive labeled risk index through risk decoding.

[0117] Specifically, step S42: Construct a defect coupling enhancement representation vector based on the multi-layer defect information map. The defect coupling enhancement representation vector is then used to obtain the comprehensive labeled risk index of each grid in the risk quantification analysis area through risk decoding.

[0118] For each geographic grid in the risk quantification analysis area, quantify the visual interference defect value. Quantitative values ​​of risk defects Combine the two elements to form an initial 2D feature vector. This serves as the initial defect description for the grid. To achieve coupling and enhancement between defects, the initial feature vectors of all grids are input into a graph neural network (GNN). Through multi-layer information propagation and aggregation mechanisms, the model can learn and capture the complex spatial dependencies between the central grid defect and its neighboring defects, thereby transforming and enhancing the initial low-dimensional vector into a high-dimensional, such as 128-dimensional, defect coupling enhancement representation vector. This vector deeply integrates comprehensive risk information about itself and its surroundings. This high-dimensional representation vector... The input consists of a risk decoding module made up of fully connected layers, which is fed through a learnable weight matrix. and bias terms A linear transformation is performed to map complex feature information into a single logical value. This logical value is then activated by the Sigmoid activation function. After normalization, a comprehensive labeled risk index between 0 and 1 is obtained. This index precisely quantifies the probability of errors in the labeled data of the corresponding grid; the closer the value is to 1, the higher the risk. This is the Sigmoid activation function.

[0119] Based on the comprehensive annotation risk index, a risk classification response and pose correction closed loop is used to make the final validity judgment and processing of the ground feature annotation;

[0120] Specifically, step S43: Based on the comprehensive annotation risk index of each grid in the risk quantification analysis area, perform risk classification response and pose correction closed loop to make the final validity judgment and processing of the ground feature annotation;

[0121] The risk-based response and pose correction closed-loop analysis includes the following:

[0122] The comprehensive annotation risk index is compared with the preset warning threshold and critical threshold to classify the ground feature annotations. If the ground feature annotation area does not belong to the risk quantification analysis area or the comprehensive annotation risk index of the risk quantification analysis area is less than the warning threshold, it is determined to be valid. If the comprehensive annotation risk index is between the warning threshold and the critical threshold, that is, the comprehensive annotation risk index is greater than or equal to the warning threshold and less than or equal to the critical threshold, a warning mark is added and manual confirmation is required. If the comprehensive annotation risk index is higher than the critical threshold, the pose smoothing extrapolation strategy is immediately triggered. The inertial measurement unit data is used to perform accurate integration calculations during periods when visual signals are unreliable to correct the UAV pose, thereby ensuring the final validity of the annotation.

[0123] Based on the comprehensive labeled risk index of each grid in the risk quantification analysis area. A risk classification response and pose correction closed loop is initiated to obtain the final feature annotation. If the feature annotation area does not belong to the risk quantification analysis area, or the comprehensive annotation risk index of the risk quantification analysis area is... <Warning Threshold This indicates a good visual positioning environment, and all labeled data are deemed qualified, becoming the final valid feature labels; if the comprehensive labeling risk index is in the warning range, that is... , If the system determines that the confidence level of the current labeled data has decreased due to a critical threshold, the label will still be displayed, but a warning label will be automatically added to prompt the operator to conduct a secondary confirmation. Only after manual confirmation can the label be considered finally valid. The risk index of the risk quantification analysis area grid... This indicates that a severe global drift has occurred or is about to occur in the Visual Odometry (VIO). The system will immediately trigger a pose smoothing extrapolation strategy, using high-frequency inertial sensor data to correct the UAV's pose, thus ensuring that ground feature annotations generated during this period can still achieve accurate positioning. These thresholds are determined based on empirical statistics, through analysis of thousands of hours of data covering both favorable and unfavorable flight environments. Favorable and unfavorable flight environments include solid-color walls and water bodies. The 95th percentile of the risk index distribution under normal operating conditions is set as... The 10th percentile of the risk index distribution under abnormal operating conditions is set as... To ensure high recall and low false alarm rate, the pose smoothing extrapolation strategy focuses on using data from the inertial measurement unit (IMU) to perform precise integration calculations during brief gaps in the unreliable visual signal to extrapolate the continuous motion trajectory of the UAV. The specific process is as follows: All unreliable visual data is paused, and the last reliable pose state before VIO failure is accurately recorded. This includes the UAV's position in the world coordinate system. ,speed and posture Attitude is typically represented by a quaternion. The system continuously acquires raw measurements from the IMU at a high frequency, namely linear accelerations along the three axes of the body coordinate system. and angular velocities along the three axes , where h represents a direct measurement derived from hardware. At each extremely short time step... Internally, the system performs integral calculations using the following kinematic equations: The acceleration measured by the IMU in the body coordinate system... Transforming to the world coordinate system involves using the pose from the previous moment. To complete this, obtain the acceleration in the world coordinate system. Because the IMU's accelerometer cannot distinguish between the drone's own acceleration and gravitational acceleration, the measured... This includes a gravitational component. To obtain the true acceleration of the drone, the gravitational acceleration vector must be subtracted. Its value is approximately [0, 0, -9.8]. Realistic motion acceleration 0. After obtaining the actual acceleration, update the velocity and position. New velocity From the velocity of the previous moment and the actual acceleration at the current moment The calculation shows that: Next, using the position from the previous moment... ,speed and current real acceleration To calculate the current new position : Simultaneously, the angular velocity measured by the IMU... Integrating data to update the drone's attitude, from... Calculate the pose at the current moment. By repeatedly performing the above integration process at a high frequency within 1 to 2 seconds, the system can generate a smooth pose trajectory that conforms to the laws of physical motion, thereby effectively filling the positioning data gap caused by VIO failure and ensuring the continuity and robustness of the system's pose output in complex environments.

[0124] Step S43 is the final execution and closed-loop stage of the entire method. It transforms all the results of the preliminary analysis and prediction into specific, actionable decisions. The solution lies in establishing a differentiated processing strategy based on risk levels. The core of this step is to use the "comprehensive annotation risk index" generated in step S42. By comparing it with preset "early warning thresholds" and "critical thresholds," risks are divided into three levels, each with distinctly different response measures. Annotations in low-risk areas are directly adopted, reflecting efficiency; annotations in medium-risk warning areas introduce a manual secondary confirmation mechanism, striking a balance between efficiency and reliability; while high-risk areas decisively trigger a "pose smoothing extrapolation strategy," suspending the adoption of unreliable visual data and instead using high-frequency IMU data for integral extrapolation to correct pose. This is the last line of defense to ensure positioning continuity and annotation accuracy in extreme cases. This step is the endpoint of the entire method's value realization and is indispensable. It constitutes a complete closed loop from risk perception to risk management. Without step S43, all the analysis, diagnosis, documentation, and prediction work done in the preceding steps (S10 to S42) would be in vain. While the system can identify where risks exist and their magnitude, it lacks any mechanism to apply this information to guide actual feature annotation. The determination of feature annotation validity remains arbitrary, unable to differentiate based on risk level, and even less capable of activating contingency plans to salvage location data when the VIO system is about to experience or has already experienced severe drift. The entire approach will merely be a complex "risk analyzer," failing to become an intelligent system capable of "real-time interpretation and annotation" and ensuring the validity of results.

[0125] Example 2

[0126] like Figure 2 As shown, a real-time interpretation and ground feature annotation system for UAV aerial survey images includes an acquisition and positioning module, a visual interference module, a spatiotemporal evolution module, and a judgment and processing module, specifically:

[0127] The data acquisition and positioning module collects and analyzes the internal data flow of multidimensional positioning solution in real time, identifies trajectory drift points using a dynamic baseline model, and determines the pose solution turbulence events and their time windows by combining clustering and correlation analysis of characteristic curves.

[0128] The visual interference module uses the determined turbulence events and their time windows to backtrack and extract the corresponding high frame rate image sequences. By analyzing the dense motion vector field and feature loss areas between images, it identifies the causal interference domain that causes the turbulence and projects it back into the physical space to generate a geospatial distribution map of the visual interference source.

[0129] The spatiotemporal evolution module determines the risk analysis area based on the geospatial distribution map of visual interference sources. It comprehensively calculates multiple risk factors within the risk analysis area, such as the lack of basic texture, changes in illumination, and the occupation of dynamic targets. Combined with the UAV's field of view footprint of the preset mission route, it generates a spatiotemporal evolution prediction map of drift risk.

[0130] The judgment and processing module performs spatiotemporal alignment and unified mapping between the geospatial distribution map of visual interference sources and the spatiotemporal evolution prediction map of drift risk, decodes the comprehensive annotation risk index of each grid, and initiates a hierarchical response mechanism to obtain the final validity judgment and processing of the ground feature annotations.

[0131] It should be understood that although the steps in the flowcharts of the various embodiments of the present invention are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the various embodiments may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least a portion of the sub-steps or stages of other steps.

[0132] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments described above. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include non-volatile and / or volatile memory. Non-volatile memory can include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM), or flash memory. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM is available in various forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), dual data rate SDRAM (DDRSDRAM), enhanced SDRAM (ESDRAM), synchronous link DRAM (SLDRAM), Rambus direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0133] The foregoing description is illustrative of the invention and should not be construed as limiting it. Although several exemplary embodiments of the invention have been described, those skilled in the art will readily understand that many modifications can be made to the exemplary embodiments without departing from the novel teachings and advantages of the invention. Therefore, all such modifications are intended to be included within the scope of the invention as defined in the claims. It should be understood that the foregoing description is illustrative of the invention and should not be construed as limiting it to the specific embodiments disclosed, and modifications to the disclosed embodiments and other embodiments are intended to be included within the scope of the appended claims. The invention is defined by the claims and their equivalents.

Claims

1. A method for real-time interpretation and feature annotation of UAV aerial survey images, characterized in that, Includes the following steps: By collecting and analyzing the internal data flow of multidimensional positioning solution in real time, the trajectory drift point is identified using the dynamic baseline model. Combined with the clustering and correlation analysis of characteristic curves, the pose solution turbulence event and its occurrence time window are determined. By utilizing the determined turbulence events and their time windows, and by analyzing the dense motion vector field and feature loss areas between images, the causal interference domain of turbulence is identified, and a geospatial distribution map of visual interference sources is generated. Based on the geospatial distribution map, the risk analysis area is determined. Taking into account the risk factors of lack of basic texture, changes in lighting, and the occupation of dynamic targets, a spatiotemporal evolution prediction map of drift risk is generated. The geospatial distribution map of visual interference sources is spatiotemporally aligned and uniformly mapped with the spatiotemporal evolution prediction map of drift risk. The comprehensive annotation risk index of each grid is decoded to obtain the final validity judgment and processing of the ground feature annotation. The internal data stream of the multidimensional positioning solution is a data stream composed of the simultaneous acquisition of the visual odometry scale factor, the determinant value of the fusion filter state covariance submatrix, and the distribution entropy of effective visual feature points, all with timestamps. The feature curves are the first visual feature point distribution entropy curve, the second filter state covariance curve, and the third scale factor curve constructed from the internal data stream of the multidimensional positioning solution. The feature loss region is the convex hull region formed by the effective visual feature points that failed to be tracked in the high frame rate image sequence before and after the occurrence of the turbulence event. The steps for determining the pose calculation turbulence event and its occurrence time window are as follows: Construct a first visual feature point distribution entropy curve, a second filtered state covariance curve, and a third scale factor curve based on the internal data flow of the multidimensional positioning calculation; obtain a set of discrete drift time points based on single-source nominal trajectory drift point analysis, and define a drift analysis window through clustering; identify the first causal instability characteristic relationship and the second symptom resonance characteristic relationship within the drift analysis window to determine the pose calculation turbulence event, and record the start time and duration of the pose calculation turbulence event. The start time is the moment when the first visual feature point distribution entropy curve shows a significant local minimum, and its duration is from the start time until the value of the second filtered state covariance curve falls back from its significant local maximum. The time elapsed until the determinant value of the fusion filter state covariance submatrix falls below the upper threshold; wherein, a data point is determined to be a significant local minimum if it meets the following two conditions: the value of the data point is lower than the values ​​of the two data points immediately preceding and following it on the time axis, forming a local low point; and the value of the data point is lower than a pre-set lower threshold for the distribution entropy of visual feature points, representing a serious deterioration in information quality; a data point is determined to be a significant local maximum if it meets the following two conditions: the value of the data point is higher than the values ​​of the two data points immediately preceding and following it on the time axis, forming a local high point; and the value of the data point exceeds a pre-set upper threshold for the determinant value of the fusion filter state covariance submatrix, representing excessively high system uncertainty; The steps of analyzing the dense motion vector field and feature loss region between images, identifying the bruising causal interference domain, and generating a geospatial distribution map of visual interference sources are as follows: using the determined pose to solve the bruising event to determine the problem time window, so as to obtain a high frame rate image sequence before and after the event occurrence time; The dense motion vector field between two consecutive frames in a high frame rate image sequence before and after the event is calculated. Motion mode decoupling and saliency discrimination are performed on the dense motion vector field to obtain the non-dominant motion manifold region. Full-spectrum filtering is performed on the images in the high frame rate image sequence before and after the event to generate a set of filtering responses. Based on the response aggregation, a texture intensity energy map is calculated. The texture intensity energy map is then subjected to low-threshold segmentation and morphological optimization to obtain the low-valley region of the full-spectrum filtering response. When the spatial overlap rate of the union of the non-dominant motion manifold region and the trough region of the full-spectrum filter response is higher than the preset causal association threshold in the feature loss region, the union is identified as the causal attribution interference domain. By obtaining the interference type and influence intensity of the causal attribution interference domain, and by back-projecting it from the image space and aggregating it into the geographic space based on the UAV pose, a geographic spatial distribution map of the visual interference source is obtained.

2. The method for real-time interpretation and feature annotation of UAV aerial survey images according to claim 1, characterized in that, The steps for obtaining the single-source nominal trajectory drift point are as follows: When the UAV performs aerial surveying missions, it collects and synchronizes the internal data stream of multidimensional positioning calculation in real time, including the visual odometry scale factor, the determinant value of the fusion filter state covariance submatrix, and the distribution entropy of effective visual feature points. Based on the internal data flow of multidimensional positioning solution, dynamic baseline models based on sliding windows are established to identify outliers that exceed a preset threshold and obtain single-source nominal trajectory drift points.

3. The method for real-time interpretation and ground feature annotation of UAV aerial survey images according to claim 1, characterized in that, The steps for identifying the first causal instability characteristic relationship and the second symptom resonance characteristic relationship are as follows: Significant local minima were detected on the distribution entropy curve of the first visual feature point, and significant local maxima were detected on the covariance curve of the second filtered state, which was determined to identify the first causal unstable feature relationship. Centered on the moment of significant local maximum detected in the first causal instability characteristic relationship, the local variance of all data points of the third scale factor curve within the time window is calculated through a fixed-length time window. If the local variance value is significantly higher than its historical average level, it is determined to be the second symptom resonance characteristic relationship.

4. The method for real-time interpretation and feature annotation of UAV aerial survey images according to claim 1, characterized in that, The steps for obtaining the spatiotemporal evolution prediction map of drift risk are as follows: A three-dimensional field-of-view cone is constructed based on a preset task file, and its intersection with the three-dimensional digital surface model of the risk quantification analysis area is calculated in three-dimensional space to obtain the UAV's field-of-view footprint.

5. The method for real-time interpretation and feature annotation of UAV aerial survey images according to claim 4, characterized in that, The steps for obtaining the spatiotemporal evolution prediction map of drift risk also include: After normalizing the illumination change risk index and the dynamic target interference risk index, they are weighted and fused with the basic texture scarcity index to calculate the total risk potential of each spatial location at different times. Based on the preset task file, the total risk potential of all locations falling within the UAV's field of view footprint is summed at each moment to obtain the instantaneous risk score at that moment. These scores are then arranged in chronological order to obtain a prediction map of the spatiotemporal evolution of drift risk.

6. The method for real-time interpretation and feature annotation of UAV aerial survey images according to claim 5, characterized in that, The steps for obtaining the dynamic target interference risk index are as follows: Risk quantification analysis areas were determined based on the geospatial distribution map of visual interference sources; The risk quantification analysis area is subjected to grid analysis, the corner density of each grid cell is calculated, and the basic texture scarcity index is obtained by normalizing and reversing the corner density.

7. The method for real-time interpretation and feature annotation of UAV aerial survey images according to claim 6, characterized in that, The steps for obtaining the dynamic target interference risk index also include: A binary shadow map is generated based on the time-series solar position and a three-dimensional model. The light change risk index is obtained by calculating the time gradient of the binary shadow map. A probabilistic dynamic target occupancy grid model is established based on historical data, and the probability and expected velocity of each grid being occupied by a target are predicted to calculate the dynamic target interference risk index.

8. The method for real-time interpretation and ground feature annotation of UAV aerial survey images according to claim 1, characterized in that, The final validity determination and processing steps for the feature annotations are as follows: By spatiotemporally aligning and uniformly mapping the geospatial distribution map of visual interference sources and the spatiotemporal evolution prediction map of drift risk, a multi-layer defect information map is obtained. Based on the multi-layer defect information map, a defect coupling enhancement representation vector is constructed. The defect coupling enhancement representation vector is then used to obtain a comprehensive labeled risk index through risk decoding. Based on the comprehensive annotation risk index, a risk classification response and pose correction closed loop is used to make the final validity judgment and processing of ground feature annotations.

9. A real-time interpretation and feature annotation system for UAV aerial survey images, characterized in that, A method for real-time interpretation and feature annotation of UAV aerial survey images according to any one of claims 1-8 includes: The data acquisition and positioning module collects and analyzes the internal data flow of multidimensional positioning solution in real time, identifies trajectory drift points using a dynamic baseline model, and determines the pose solution turbulence events and their time windows by combining clustering and correlation analysis of feature curves. The visual interference module utilizes the determined turbulence events and their time windows to identify the causal interference domain of turbulence by analyzing the dense motion vector field and feature loss areas between images, and generates a geospatial distribution map of visual interference sources. The spatiotemporal evolution module determines the risk analysis area based on the geospatial distribution map, and generates a spatiotemporal evolution prediction map of drift risk by taking into account risk factors such as lack of basic texture, changes in lighting, and the occupation of dynamic targets. The judgment and processing module performs spatiotemporal alignment and unified mapping between the geospatial distribution map of visual interference sources and the spatiotemporal evolution prediction map of drift risk, and decodes the comprehensive annotation risk index of each grid to obtain the final validity judgment and processing of the ground feature annotation.

Images