Space-air-ground coordinated natural resource three-dimensional investigation and monitoring method, device, equipment and medium

Abstract

The application relates to a space-air-ground coordinated natural resource three-dimensional investigation and monitoring method, device, equipment and medium. The method converts a digital elevation model into a navigation topological manifold to utilize terrain constraints as navigation priors, establishes a space-air-ground three-element coupling optimization model in combination with satellite orbit prediction, and generates unmanned aerial vehicle (UAV) cooperative planning by adopting distributed predictive control, so that a UAV group is autonomously reconstructed into an optimal observation array based on terrain topology within a satellite revisit window, positioning and formation keeping under a global navigation satellite system (GNSS) signal are realized through a terrain fingerprint atlas, and finally, through space-air-ground data fusion and multi-period change detection, high-precision and high-efficiency cooperative monitoring of natural resource investigation under complex terrain is realized.

Description

Technical Field

[0001] This invention belongs to the field of integrated air-space-ground collaborative monitoring technology, and in particular relates to a method, device, equipment and medium for integrated air-space-ground three-dimensional survey and monitoring of natural resources. Background Technology

[0002] Currently, the field of natural resource surveys and monitoring mainly relies on a technical system combining satellite remote sensing, aerial photography, and ground surveys. However, it faces significant technical bottlenecks in complex terrain areas such as canyons and sinkholes. Traditional UAV swarm operations typically treat terrain as a navigation obstacle, employing obstacle avoidance strategies that lead to frequent maneuvers. This not only increases energy consumption but also causes formation instability and disrupts mission timing. Especially when satellite revisit windows are extremely limited, this passive obstacle avoidance mode struggles to complete full-domain data collection within the specified time.

[0003] Existing research on GNSS denied environments largely focuses on single-unit visual / inertial navigation fusion positioning, failing to fully explore the navigational value of terrain features. Although some studies have attempted to use topographic maps to assist positioning, a dynamic correlation mechanism between terrain features and satellite imaging geometry has not yet been established. When satellites pass overhead, UAV swarms often fail to form an optimal observation formation due to a lack of coordinated planning, resulting in spatial misalignment between acquired data and satellite imagery, severely impacting the accuracy of subsequent data fusion.

[0004] In terms of swarm collaboration, existing methods often separate obstacle avoidance from formation control, lacking unified spatial topology modeling. Frequent starts, stops, and turns by UAVs in complex terrain not only shorten effective operating time but also exacerbate the burden on inter-UAV communication. Furthermore, traditional multi-sensor fusion schemes fail to transform terrain occlusion effects into usable navigation constraints, losing GNSS signals in terrains such as canyons and failing to fully utilize natural geographical guidance features.

[0005] Collaborative operations between satellites and drones are still in their early stages, with each often planning its mission independently. The fixed orbital period of satellites and the dynamic operational rhythm of drones inherently conflict. When the satellite's side-swing angle changes, the ground-based drone swarm cannot adjust its observation attitude in real time, resulting in systematic geometric discrepancies between air-space and ground-based data. This spatiotemporal asynchrony is particularly prominent in natural resource monitoring that requires multi-period comparisons.

[0006] Therefore, there is an urgent need for a new survey method that can transform terrain constraints into navigation advantages and achieve dynamic coordination between air, land, and space, in order to break through the technical barriers in complex environments and improve the overall effectiveness of three-dimensional natural resource surveys. Summary of the Invention

[0007] Therefore, it is necessary to provide a method, device, equipment, and medium for three-dimensional survey and monitoring of natural resources that integrates air, space, and ground to address the aforementioned technical problems.

[0008] Firstly, this application provides a three-dimensional survey and monitoring method for natural resources that integrates air, space, and ground, including:

[0009] S1. Construct a navigation topology manifold based on the digital elevation model data of the target area; wherein, the navigation topology manifold is used to transform terrain constraints into prior information for UAV swarm navigation;

[0010] S2. Based on satellite orbital parameters and satellite imaging geometry, predict the satellite's revisit time window and ground coverage area using the Kepler orbital equation; based on the revisit time window and ground coverage area, establish a ternary coupled optimization model of satellite-terrain-UAV swarm in conjunction with the navigation topology manifold;

[0011] S3. Solve the ternary coupled optimization model through a distributed model predictive control framework to generate a collaborative task planning scheme for the UAV swarm; the collaborative task planning scheme includes the waypoint sequence, velocity profile and task time window for each UAV.

[0012] S4. Within the revisit time window, control the UAV swarm to enter the target area based on the collaborative mission planning scheme; perform topological self-reconstruction based on the navigation topology manifold to form an observation array that is geometrically matched with satellite imaging;

[0013] S5. Extract terrain feature points in the target area and construct a terrain fingerprint map based on the terrain feature points; based on the terrain fingerprint map, control the UAV swarm to perform relative navigation and formation maintenance in a GNSS denied environment through collaborative simultaneous localization and map building.

[0014] S6. At the moment when the satellite passes overhead, the UAV swarm controlled by the observation array synchronously collects surface data, and fuses the surface data with the satellite imagery to obtain air-space-ground integrated data; based on the change detection results and trend prediction results of the multi-period air-space-ground integrated data, the results of the three-dimensional survey and monitoring of natural resources are generated.

[0015] Secondly, this application also provides a three-dimensional survey and monitoring device for natural resources that integrates air, space, and ground, for implementing the method described in the first aspect. The device includes:

[0016] The terrain navigation manifold modeling module is used to construct a navigation topology manifold based on the digital elevation model data of the target area; the navigation topology manifold is used to transform terrain constraints into prior information for UAV swarm navigation.

[0017] The satellite-terrain co-optimization module is used to predict the satellite's revisit time window and ground coverage area based on satellite orbital parameters and satellite imaging geometry using the Kepler orbital equations; and to establish a ternary coupled optimization model of satellite-terrain-UAV swarm based on the revisit time window and ground coverage area, combined with the navigation topology manifold.

[0018] The collaborative task planning module solves the ternary coupled optimization model through a distributed model predictive control framework to generate a collaborative task planning scheme for the UAV swarm. The collaborative task planning scheme includes the waypoint sequence, velocity profile and task time window for each UAV.

[0019] The UAV formation self-reconstruction module is used to control the UAV swarm to enter the target area based on the collaborative mission planning scheme within the revisit time window; it performs topological self-reconstruction based on the navigation topology manifold to form an observation formation that is geometrically matched with satellite imaging.

[0020] The terrain feature extraction and navigation module is used to extract terrain feature points in the target area and construct a terrain fingerprint map based on the terrain feature points. Based on the terrain fingerprint map, through collaborative simultaneous localization and map building, the module controls the UAV swarm to perform relative navigation and formation maintenance in GNSS-denied environments.

[0021] The air-space-ground data fusion module is used to control a swarm of drones to synchronously collect surface data based on the observation array when the satellite passes overhead. The surface data is then fused with satellite imagery to obtain air-space-ground fused data. Based on the change detection results and trend prediction results of the multi-period air-space-ground fused data, the results of the three-dimensional survey and monitoring of natural resources are generated.

[0022] Thirdly, this application also provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement a three-dimensional survey and monitoring method for natural resources that integrates air, space, and ground as described in the first aspect.

[0023] Fourthly, this application also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements a three-dimensional survey and monitoring method for natural resources that combines air, space, and ground operations as described in the first aspect.

[0024] The aforementioned method, device, equipment, and medium for three-dimensional survey and monitoring of natural resources using a combined air-space-ground approach transforms digital elevation models into topological manifolds with navigational guidance value. This transforms terrain constraints, traditionally considered obstacles, into favorable conditions for UAV path planning. Based on satellite orbital mechanics prediction revisit windows, a coupled optimization model integrating air, space, and ground is established. Distributed model predictive control enables collaborative task planning for UAV swarms. Within the limited time window of satellite overpasses, high-precision fingerprint maps are constructed using terrain features. Multi-UAV collaborative SLAM technology achieves relative positioning and stable formation maintenance in complex terrain environments where GNSS signals are completely rejected. Finally, through deep fusion of multi-source air-space-ground data and temporal change detection and analysis, the method effectively solves the technical bottlenecks of traditional survey methods in complex terrain environments, such as poor positioning accuracy, low data acquisition efficiency, and difficulties in air-space-ground collaboration. This represents a technological leap in natural resource survey and monitoring, from single-method to three-dimensional collaboration, from passive obstacle avoidance to active terrain utilization, and from decentralized operations to system optimization. It significantly improves the quality of data acquisition and monitoring efficiency in complex environments. Attached Figure Description

[0025] To more clearly illustrate the technical solutions in the embodiments or related technologies of this application, the accompanying drawings used in the description of the embodiments or related technologies will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0026] Figure 1 A flowchart illustrating a three-dimensional survey and monitoring method for natural resources that integrates air, space, and ground operations, provided by this invention.

[0027] Figure 2 This is a schematic diagram illustrating the process of constructing a navigation topology manifold in one optional embodiment of the present invention;

[0028] Figure 3 This is a schematic diagram of the structure of a three-dimensional natural resource survey and monitoring device that integrates air, space, and ground, provided by the present invention. Detailed Implementation

[0029] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0030] refer to Figure 1 The application presents a flowchart illustrating a three-dimensional survey and monitoring method for natural resources that integrates air, space, and ground, comprising the following steps:

[0031] S1. Construct a navigation topology manifold based on the digital elevation model data of the target area; wherein, the navigation topology manifold is used to transform terrain constraints into prior information for UAV swarm navigation.

[0032] Specifically, the acquisition of Digital Elevation Model (DEM) data for the target area is based on the accuracy requirements of the survey and monitoring, as well as the characteristics of the regional scope, selecting an appropriate data source. For high-precision monitoring scenarios at the micro- and meso-scale, DEM data acquired through aerial photogrammetry or UAV oblique photogrammetry can be selected; for macro-scale, large-scale monitoring tasks, DEM data generated by satellite remote sensing technology can be used. After acquiring the raw DEM data, preprocessing operations are performed, mainly including data denoising and data normalization. The data denoising stage employs a filtering algorithm adapted to the DEM data resolution to eliminate random noise introduced by sensor errors, external interference, and other factors in the raw data, ensuring the accuracy of the terrain elevation information. A Gaussian filtering algorithm can be used here, with the following formula: ,in Indicates the filter kernel in coordinates The weight value at that point, The standard deviation of the Gaussian function is used to control the smoothness of the filter. and These represent the horizontal and vertical coordinates of pixels within the filter kernel relative to the kernel center, respectively; data normalization transforms DEM data from different data sources and coordinate systems into a preset geographic coordinate framework, using the following normalization formula: ,in This indicates the normalized DEM at the pixel level. The elevation value at that location, Indicates the original DEM in pixels The elevation value at that location, This represents the minimum elevation value in the original DEM data. This represents the maximum elevation value in the original DEM data. This formula can be used to map the elevation value to the [0,1] interval, laying the foundation for subsequent terrain feature analysis and topology construction.

[0033] After DEM data preprocessing, the next stage is terrain feature extraction. Using terrain feature recognition algorithms, key features reflecting terrain undulation and structure are extracted from the DEM data, including but not limited to ridgelines, valley lines, abrupt slope changes, and significant elevation control points. Here, a slope calculation algorithm is used to identify areas of abrupt slope changes. The slope calculation formula is as follows: ,in Indicates DEM in pixels The slope value at that location, , Representing pixels Elevation values ​​of adjacent pixels, , Representing pixels Elevation values ​​of the left and right adjacent pixels, and These represent the horizontal and vertical resolutions of the DEM data, i.e., the actual distance between adjacent pixels on the ground. These terrain features are core elements for constructing the navigation topology manifold and can accurately characterize the terrain constraints of the target area.

[0034] The process of constructing the navigation topology manifold then begins. First, the basic building blocks of the topology manifold are defined: topology nodes and topology edges. Topology nodes are selected from extracted terrain feature points, and each topology node contains not only its spatial location information in the geographic coordinate system. ( , They represent the first The planar coordinates of the topological nodes, This indicates its elevation coordinates and also integrates terrain attribute parameters such as elevation, slope, and aspect corresponding to that location. Topological edges are used to connect adjacent topological nodes, and their construction must meet two core conditions: first, the terrain profile between adjacent topological nodes must meet the safety constraints of UAV flight, i.e., there must be no angle of ascent exceeding the UAV's maximum. Maximum descent angle Restricted terrain obstacles, adjacent nodes and The formula for calculating the climb angle between them is: The formula for calculating the descent angle is: , must meet and Secondly, the distance between adjacent topology nodes must be within the effective coverage area of ​​the UAV swarm communication link. Inside, adjacent nodes and The formula for calculating the distance between them is: , must meet This ensures that the drone can maintain stable information exchange while flying along the topological edge.

[0035] After defining the topological nodes and edges, a complete navigation topological manifold is constructed using a topology generation algorithm. This topological manifold is essentially a spatial topological structure that reflects terrain constraints and the feasibility of UAV flight, transforming complex terrain information within the target area into prior information usable for UAV swarm navigation. During subsequent UAV swarm navigation, the navigation system can quickly eliminate infeasible flight paths caused by terrain obstacles based on this topological manifold, significantly reducing the computational complexity of navigation path planning while ensuring terrain safety during UAV swarm flight and avoiding flight accidents caused by terrain occlusion or obstacle collisions.

[0036] S2. Based on satellite orbital parameters and satellite imaging geometry, predict the satellite's revisit time window and ground coverage area using the Kepler orbital equation; based on the revisit time window and ground coverage area, establish a ternary coupled optimization model of satellite-terrain-UAV swarm in conjunction with the navigation topology manifold.

[0037] Specifically, satellite orbital parameters mainly include the semi-major axis of the orbit. Orbital eccentricity Track inclination Right ascension of ascending node Perigeal argument and the near point angle These parameters can be obtained from orbital data files released by the satellite mission or orbital prediction data provided by the aerospace telemetry and control department, and can completely describe the satellite's trajectory in the celestial coordinate system. Satellite imaging geometric parameters cover the field of view of the satellite imaging sensor. Imaging resolution Side swing angle Range and optical axis orientation at imaging time (Unit vector), these parameters are determined by the satellite's design specifications and mission planning, and directly determine the satellite's imaging coverage capability and geometric projection relationship of ground target areas.

[0038] Based on the acquired satellite orbital parameters, the Kepler orbital equations are used to predict the satellite's revisit time window and ground coverage area. The expression for the Kepler orbital equations is as follows: ,in The mean anomaly angle, physically speaking, is the angle a satellite rotates through as it moves from its perigee to its current position with an average angular velocity. It can be found in the initial mean anomaly angle in the satellite's orbital parameters. With the average angular velocity of the orbit (Depend on calculate, (the Earth's gravitational constant), combined with the time variable. (The time interval from the initial time to the current time) is expressed by the formula Calculated; The deviated anomaly angle is an intermediate variable introduced to simplify the solution of Kepler's equations, used to establish the relationship between the mean anomaly angle and the true anomaly angle of the satellite. The relationship between them, true nearest angle Angle with near point The conversion relationship is as follows ; Indicates orbital eccentricity, used to describe the degree of ellipticity of a satellite's orbit. The orbit is circular. and The orbit is elliptical.

[0039] By solving the Kepler orbital equations using numerical iterative algorithms (such as the Newton-Raphson iterative method), the true anomaly angle of the satellite at any given time can be obtained. The satellite's spatial position vector in the celestial coordinate system is calculated by combining the satellite's orbital parameters. With velocity vector The formula for calculating the position vector is: ,in and These are unit orthogonal vectors in the orbital plane. Subsequently, a coordinate transformation algorithm (converting sequentially from the celestial coordinate system to the geocentric Earth-fixed coordinate system, and then from the geocentric Earth-fixed coordinate system to the local coordinate system of the target area) is used to transform the satellite's celestial coordinate system to the Earth-fixed coordinate system. Then, combined with the satellite's imaging geometric parameters, the coverage area of ​​the ground by the satellite at different times is calculated. The boundary of the satellite's ground coverage area can be calculated using the field of view of the imaging sensor, and the radius of the ground coverage circle is... The calculation formula is ,in The orbital altitude of the satellite above the target area. The field of view of the imaging sensor is represented by a polygon with latitude and longitude coordinates, which is the satellite's ground coverage area. Simultaneously, based on the periodicity and orbital return characteristics of the satellite's orbit around the Earth, the pattern of repeated coverage of the target area by the satellite is analyzed to determine the time intervals within which the satellite will pass over the target area multiple times in the future. ( (This refers to the number of revisits), and this time interval is the satellite's revisit time window. Each revisit time window corresponds to one effective imaging opportunity for the satellite of the target area.

[0040] After obtaining the satellite revisit time window and ground coverage area, a three-element coupled optimization model of satellite-terrain-UAV swarm is established based on the navigation topology manifold constructed in step S1. The objective function of this optimization model is set to maximize the effective data acquisition area and data acquisition accuracy of the UAV swarm within the satellite revisit time window, while minimizing the energy consumption and mission completion time of the UAV swarm. The mathematical expression of the objective function is as follows: ,in The objective function value, , , , The weighting coefficients for effective data acquisition area, data acquisition accuracy, total energy consumption, and total task time are respectively (satisfying...) ), The effective data collection area of ​​the UAV swarm (i.e., the collection area within the satellite ground coverage area and meeting terrain constraints). The accuracy of data acquisition is determined by the degree of matching between the UAV swarm observation array and the satellite imaging geometry; the higher the degree of matching, the greater the accuracy value. The total energy consumption of the drone swarm (calculated from the flight path length and flight speed of each drone), , For the number of drones, For the first A drone at all times With speed Power consumption during flight , The first (Start and end times of the drone mission) The total mission completion time of the drone swarm (i.e., the time difference between the time the last drone completes its mission and the time the first drone starts its mission).

[0041] The model's constraints fall into three categories: The first category is terrain constraints, directly provided by the navigation topology manifold, which restrict the flight path of the UAV swarm to within the feasible region defined by the topology manifold, i.e., the first... Location of the drone Must meet ( The first type is the feasible region of the navigation topology manifold, which requires that the flight meet constraints such as terrain elevation and slope to avoid terrain collisions; the second type is satellite constraints, which limit the observation area of ​​the UAV swarm to within the satellite ground coverage area. , (For satellite ground coverage), data acquisition time must be within the satellite revisit time window. , For the first The data acquisition time of each drone is considered, and the geometric parameters of the drone swarm's observation formation must match the projection geometry parameters of the satellite image to ensure spatial alignment between the acquired data and the satellite image; the third category is the constraints of the drone swarm itself, including the maximum flight speed of each drone. ( Maximum battery life ( Maximum climb / descent angle , (The calculation method for the ascent and descent angles is consistent with that in S1, and must meet the following requirements) and Performance constraints such as sensor operating range, and communication distance constraints within the drone swarm (distance between adjacent drones). ), formation safety distance constraints ( , Measures such as safe distance thresholds are used to ensure that drone swarms can stably and safely perform collaborative tasks.

[0042] S3. Solve the ternary coupled optimization model through a distributed model predictive control framework to generate a collaborative task planning scheme for the UAV swarm; the collaborative task planning scheme includes the waypoint sequence, velocity profile and task time window for each UAV.

[0043] Specifically, the first step is to build a distributed model predictive control (DMPC) framework. The core advantage of this framework is that it can decompose the complex UAV swarm cooperative optimization problem into multiple sub-problems, which are solved in a distributed manner by the local controller corresponding to each UAV. At the same time, global coordination is achieved through information interaction between UAV swarms, which effectively reduces the computational complexity under the centralized control framework and improves real-time performance and fault tolerance.

[0044] The distributed model predictive control framework comprises three key layers: First, the global coordination layer, responsible for receiving the objective function and constraints of the ternary coupled optimization model, formulating global collaborative optimization strategies, sending global optimization indicators and collaborative constraint information to each UAV's local controller, and simultaneously receiving local optimization results from each local controller to evaluate and adjust the global optimization effect; second, the local control layer, where each UAV is equipped with an independent local controller, based on information from the global coordination layer and its own state information (such as current position). ,speed Remaining battery power (etc.), decomposing the ternary coupled optimization model into local optimization subproblems for a single UAV, the first... The local optimization objective function of the UAV is ,in , , , The first The system includes three layers: the effective data acquisition area, data acquisition accuracy, energy consumption, and mission time of the UAVs; the third layer is the information interaction layer, which realizes information interaction between the global coordination layer and the local control layer, and between the local control layers, through the self-organizing communication network of the UAV swarm. The transmitted information includes UAV status data, local optimization results, and cooperative constraint parameters. The Time Division Multiple Access (TDMA) communication protocol is used to avoid information conflicts and ensure the real-time performance and reliability of communication.

[0045] Within the distributed model predictive control framework, a rolling time-domain optimization strategy is employed to solve the ternary coupled optimization model. This strategy divides the entire mission cycle of the UAV swarm into multiple consecutive finite time domains (i.e., prediction time domains). Within each prediction time domain, each local controller bases its prediction on the current time. System status ( For the first The system evolution model is constructed by predicting the system's evolution process in the future prediction time domain using the state vector of the UAV. The model uses the UAV's position, velocity, and remaining battery power as core state variables, and its continuous-time dynamic equation is:

[0046]

[0047]

[0048]

[0049] in, For the first The rate of change of position (i.e., velocity vector) of the drone. The unit vector for the flight direction (its value ranges to a unit sphere in three-dimensional space); The rate of change of velocity (i.e., acceleration) (This requires meeting the constraints of the UAV power system) ( For the first (Maximum acceleration of the drone). The rate of change of remaining electricity. For power consumption functions, empirical models are typically used. ( This refers to standby power. , The power coefficient, which is related to speed, is calibrated from the parameters of the UAV's power system.

[0050] To adapt to numerical computation, the continuous-time model is discretized, and a discrete time step is set. (The time interval can be taken as 0.1-1 seconds, adjusted according to the real-time requirements of the task), then the discretized state transition equation is:

[0051]

[0052]

[0053]

[0054] in, ( , To predict the start time, (for predicting the number of time-domain steps) , , The first The predicted location, speed, and remaining battery power of the step.

[0055] Based on the above system evolution model, each local controller uses the interior-point method or sequential quadratic programming method to solve the local optimization subproblems. In addition to the terrain, satellite, and self-performance constraints mentioned above, discretized state constraints also need to be added: ( (To achieve the minimum flight speed and avoid stalling) ( (Set the minimum remaining battery power threshold to ensure safe return).

[0056] In the waypoint sequence generation stage, candidate waypoints are first screened from the node set of the navigation topology manifold using the A* heuristic search algorithm. The heuristic function is set as follows:

[0057]

[0058] in, For the current topology node, For the target observation node, Let be the Euclidean distance between nodes. The terrain energy consumption coefficient at the node (determined by slope and elevation difference). , These are the weighting coefficients. The resulting candidate waypoint sequences are selected. ( (The number of waypoints) is further adjusted through nonlinear optimization, with the optimization objective being to minimize the total energy consumption and path smoothness of the path.

[0059]

[0060] in, For the first Cruising speed of the route The first term is the smoothness weight, and the second term is the path curvature penalty term, which ensures the drone's flight attitude is stable.

[0061] Velocity profile generation is based on a waypoint sequence, employing a dynamic programming algorithm to divide flight phases (acceleration, cruise, deceleration). Each flight segment... length Velocity is the distance dimension of the state space. As a velocity dimension, state is defined. The cost is from the starting point to the... The section of the route is based on speed The cumulative energy consumption of flight is expressed by the state transition equation:

[0062]

[0063] in, For the first The average speed of the route Indicates the first The length of the route For the first The section of the route is based on speed The cumulative cost of flight. By backtracking, the optimal speed for each segment of the route can be obtained, thus forming a complete speed profile including acceleration time, cruising speed, and deceleration time.

[0064] Determining the mission time window requires combining the satellite revisit time window with flight time predictions, and applying this to each waypoint. Its earliest arrival time With the latest arrival time They are respectively:

[0065]

[0066]

[0067] in, For the first Task completion time for each waypoint , The first Maximum and minimum permissible speeds for the flight path segment. Combined with satellite revisit time windows. Constraints were adjusted on the time windows of all waypoints to ensure that the time windows of the waypoints corresponding to the data acquisition tasks fell entirely within the satellite revisit time window. ( (This refers to the data acquisition time for the j-th waypoint).

[0068] The latest arrival time needs to take into account both the satellite revisit window constraint and the redundancy of subsequent tasks. The formula is:

[0069]

[0070] in, This is the end time of the satellite revisit window. To ensure all subsequent routes are at minimum speed Total flight time (to ensure subsequent missions can be completed within the window). This is the longest time that the current flight path can be traveled at the minimum speed. Therefore, the... The task time window for each waypoint is defined as If the time windows of adjacent waypoints overlap ( If this occurs, the local controller will re-optimize the speed profile, adjusting the cruise speed balance time allocation of the intermediate route to ultimately form a conflict-free mission time window sequence.

[0071] Through the above process, the distributed model predictive control framework outputs the waypoint sequence (including 3D coordinates and terrain attributes), velocity profile (including velocity-time curves for acceleration / cruise / deceleration phases), and mission time window (including time constraints for each waypoint) for each UAV, which together constitute a complete collaborative mission planning scheme.

[0072] S4. Within the revisit time window, control the UAV swarm to enter the target area based on the collaborative mission planning scheme; perform topological self-reconstruction based on the navigation topology manifold to form an observation array that is geometrically matched with satellite imaging.

[0073] Specifically, the control of the drone swarm entering the target area adopts a "position-velocity dual closed-loop feedback" mechanism. First, the flight controller of each drone collects its own position in real time. With speed Calculate the deviation from the expected trajectory:

[0074]

[0075] in, , In the collaborative task planning scheme The desired position and velocity at each moment (obtained by interpolation of the waypoint sequence and velocity profile). A PID control law is designed based on the deviation.

[0076]

[0077] In the formula, It is a proportional coefficient (adjusting the convergence speed of position deviation). This is the integral coefficient (to eliminate steady-state error). The differential coefficient (to suppress overshoot). This is for control output (used to adjust the UAV's thrust and attitude angle). During control, terrain constraints based on the navigation topology manifold are considered when the UAV approaches the topological boundary ( and distance , When the boundary threshold is used, a terrain exclusion term is introduced. ( (This is the normal vector of the topological boundary, pointing inside the feasible region) to prevent the drone from deviating from the safe flight area.

[0078] The triggering conditions for topology self-reconfiguration are divided into two categories: one is when the position deviation exceeds the limit ( , The position deviation threshold is determined by the satellite imaging resolution and can be taken as 1 / 2 of the ground distance corresponding to the imaging pixel; the second is the terrain topology change (the rate of change of local node density of the navigation topology manifold). , for Local topological region at time t, (This refers to the rate of change threshold). Once triggered, the self-reconfiguration process consists of three steps:

[0079] The first step is array target optimization. The objective function is the projection matching degree of the satellite imaging geometry.

[0080]

[0081] in, For the number of drones, The projected position of the UAV on the satellite imaging plane (pointed from the satellite's optical axis). 3D position calculation with UAVs: , (The projection matrix from the satellite imaging coordinate system to the plane coordinate system). The desired projection position (based on the satellite field of view). Uniform division: , , For the width and height of the imaging plane, (Number of rows and columns in the formation) This is the length of the diagonal of the imaging plane. The value range is [0,1]. The closer the value is to 1, the higher the matching degree between the array pattern and the satellite imaging geometry.

[0082] The second step is topology path replanning. This involves resetting the feasible node set based on the navigation topology manifold. The Hungarian algorithm is used to assign target formation positions to each UAV. And search for the optimal conversion path:

[0083]

[0084] in, For adjacent topology nodes ( and The distance, The obstacle penalty coefficient, For nodes With obstacle area distance ( Indicates accessibility. (Penalty triggered at the time).

[0085] The third step is the synchronized formation transition. The drone swarm exchanges position and speed information through distributed communication, employing a "leader-follower" cooperative control system. The drone at the center of the formation is selected as the lead drone (its trajectory is...). ), Follow the machine The expected position is ( (Assuming the relative fixed positions of the follower and navigator aircraft), by adjusting the control law Achieve synchronous conversion until the positional deviation of all drones is correct. ( (Assuming an array accuracy threshold), the observation array is constructed to geometrically match the satellite imaging.

[0086] S5. Extract terrain feature points in the target area and construct a terrain fingerprint map based on the terrain feature points; based on the terrain fingerprint map, control the UAV swarm to perform relative navigation and formation maintenance in a GNSS-denied environment through coordinated simultaneous localization and map building.

[0087] Specifically, terrain feature point extraction requires combining multi-sensor data (LiDAR point cloud, optical imagery) from the UAV, prioritizing feature points with both geometric saliency and stability, and employing an improved ORB (Oriented FAST and Rotated BRIEF) algorithm for real-time extraction. This algorithm consists of two steps: corner detection and descriptor generation. The corner detection stage uses the FAST algorithm, comparing pixel points... The corner point is determined by the difference in grayscale values ​​between the corner point and its 16 neighboring pixels. The determination formula is as follows: (in For pixels grayscale value, For neighboring pixels grayscale value, (This is a grayscale difference threshold, dynamically adjusted based on image contrast). If the number of neighboring pixels meeting this condition exceeds 8, then it is determined that... Candidate corner points are selected; then, non-maximum suppression is used to eliminate redundancy between adjacent candidate corner points, and the corner point with the largest response value is retained as the final terrain feature point. The response value is calculated using the Harris corner response function. (in The pixel grayscale covariance matrix, , , They are respectively , Gray-scale gradient in direction, For a 3×3 neighborhood window, (Empirical coefficient).

[0088] After feature point extraction, a terrain fingerprint map is constructed. This map uses "feature point spatial coordinates - local descriptor - terrain attribute" as its core structure and employs a hash table for storage and indexing. For each terrain feature point... Its atlas entries are defined as follows: ,in The three-dimensional coordinates of the feature points in the local coordinate system (obtained by direct measurement with LiDAR or stereo matching of optical images). For ORB descriptors, a 256-bit binary vector is generated by comparing the grayscale of an 8×8 pixel block in the neighborhood of the feature point, which is used for feature matching. The terrain attribute (slope) corresponding to the feature point Slope aspect ),slope The slope aspect is obtained using the slope calculation formula in S1. Through formula Calculations are performed (values ​​range from 0-360°, with 0° corresponding to true north). To improve map matching efficiency, a KD-tree (K-Dimensional Tree) is used to index the feature point coordinates. During the construction process, [the following is used as a reference]. , , The space is divided sequentially by coordinates to form a multi-level index structure. During matching, the feature points to be matched are calculated using Euclidean distance. Feature points in the map The similarity is calculated using the formula: (in for and Hamming distance, The value ranges from 0 to 1, with values ​​closer to 1 indicating a higher degree of matching.

[0089] Co-SLAM based on terrain fingerprint maps achieves global mapping through distributed information exchange by a swarm of drones. Figure 1 Consistency and individual positioning accuracy. First, establish a multi-machine information exchange protocol, where each UAV (denoted as...)... ) will use the feature points and pose estimates it observes ( The information is packaged into packets containing 3D position and 3D attitude angles, along with map fragments. These packets are then broadcast at a specific frequency to neighboring drones (drones within the communication radius) via a mesh self-organizing network. Simultaneously, the network receives information packets from neighboring drones. The information fusion stage employs a graph optimization-based approach, treating each drone's pose and feature points as nodes in a graph, and observation constraints as edges, to construct a global optimization graph. The objective function is to minimize the reprojection error of the pose and feature points.

[0090]

[0091] in, For the number of drones, for The observed set of feature points, For projection functions (projecting 3D feature points) pass position Projected onto a two-dimensional image plane. for exist Pixel coordinates of the image plane Let Euclidean distance be the square of the distance. The optimization problem is solved iteratively using the Gauss-Newton method to update the pose estimates of each UAV and the global map, thus achieving relative navigation in GNSS-denied environments.

[0092] Formation-keeping control is based on a distributed PID controller designed for relative pose error, for each UAV. Select its neighboring drones Using the formation partners determined by the observation array as a reference, calculate the relative pose error. (in for and The expected relative pose is decomposed into position error. With attitude error The position control law is:

[0093]

[0094] The attitude control law is:

[0095]

[0096] in for The expected linear velocity, For the desired angular velocity, , , These are the proportional, integral, and derivative coefficients of the position PID controller. , , These are the coefficients of the attitude PID (calibrated based on UAV dynamic parameters, typically optimized using particle swarm optimization). They are adjusted in real-time. and This ensures that the drone swarm maintains a stable observation formation in GNSS-denied environments.

[0097] S6. At the moment when the satellite passes overhead, the UAV swarm controlled by the observation array synchronously collects surface data, and fuses the surface data with the satellite imagery to obtain air-space-ground integrated data; based on the change detection results and trend prediction results of the multi-period air-space-ground integrated data, the results of the three-dimensional survey and monitoring of natural resources are generated.

[0098] Specifically, the synchronous acquisition of surface data at the moment of satellite overpass needs to achieve synchronization in both time and space. Time synchronization adopts BeiDou / GPS dual-mode time calibration. For example, the ground control station sends a time synchronization command to the drone swarm 10 minutes before the satellite overpass, calibrating the local clock of each drone to UTC time. Synchronization error control is achieved through a timestamp feedback mechanism. Each drone sends a timestamped status message to the control station, and the control station calculates the message transmission delay. ( For the time the drone sends, The control station receives the data (time) and issues compensation commands to the drones, ultimately ensuring that the time synchronization error of the drone swarm is less than or equal to a preset value. Spatial synchronization is based on the position constraints of the observation array. One minute before the satellite passes overhead, the control station sends a position fine-tuning command to the drone swarm. Through the formation-keeping control law in S5, the deviation between the actual position and the desired observation position of each drone is less than or equal to a preset value, ensuring that the collected surface data (LiDAR point cloud, hyperspectral imagery) completely overlaps with the satellite imaging coverage area in space.

[0099] The fusion of surface data and satellite imagery employs a two-step process: "feature-level registration - pixel-level fusion." Feature-level registration uses terrain feature points extracted by S5 as a bridge. First, corresponding feature points are extracted from the satellite imagery (using the ORB algorithm consistent with the UAV imagery). Then, mismatched points are eliminated using the RANSAC (Random Sample Consensus) algorithm. The number of RANSAC iterations is [not specified]. Through formula calculate( Confidence level; The ratio of the outer point; (where is the number of feature points sampled each time), and matching pairs with more than 80% of the total matching points can be retained. Based on the interior point matching pairs, the spatial transformation matrix between satellite imagery and UAV surface data is solved. (Affine transformation matrix, applicable to small areas). Solve using the least squares method The elements make satellite image feature points With drone feature points satisfy , Pixel-level fusion employs differentiated strategies for different data types. For LiDAR point clouds and satellite DSM (Digital Surface Model), a weighted average fusion is used, and the elevation values ​​of the fused DSM are... ( , (Weights are determined based on data accuracy); For hyperspectral and satellite multispectral images, wavelet transform-based fusion is used to superimpose the high-frequency details (texture, edges) of the UAV hyperspectral image with the low-frequency contours (large coverage) of the satellite multispectral image to generate a fused image with both high spatial and hyperspectral resolution.

[0100] In the change detection of multi-phase air-space-ground integrated data, two phases of integrated data are first selected. The elevation difference method is used to calculate the change area of ​​the DSM data, and the change detection threshold is set. The Otsu algorithm is used to adaptively determine (maximizing the inter-class variance between foreground and background), when... At that time, it was determined to be an area of ​​elevation change ( (This represents the elevation value of the DSM after fusion at time t); spectral differences are calculated using the spectral angle mapper (SAM) on the image data, and the SAM angle... ( , (where is the spectral vector of the two images), when ( When the angle is 10°-15°, it is determined to be the spectral variation region.

[0101] Trend forecasting is based on time series analysis and uses the ARIMA (Autoregressive Integrated Moving Average) model. The model form is as follows: ,in Let the order be the autoregressive order. It is the difference order. This refers to the moving average order. First, it involves analyzing multi-period area change data (such as monthly vegetation cover change area). Perform a stationarity test (ADF test). If the data is non-stationary, then... The difference is performed to ensure the stationarity of the differencing data; this is determined using the ACF (autocorrelation function) and PACF (partial autocorrelation function). , Value (ACF truncation corresponding) PACF truncation corresponds to ); Construct an ARIMA model and solve for the model parameters using maximum likelihood estimation. The model prediction formula is as follows: ,in This represents the predicted area. These are the autoregressive coefficients. The moving average coefficient is... For the error term, To predict the number of steps.

[0102] The final results of the three-dimensional natural resource survey and monitoring include three core categories: first, spatial results, including natural resource type distribution maps (generated based on supervised classification of fused imagery) and thematic maps of change areas (labeling change types such as vegetation degradation and water body expansion); second, attribute results, including statistical reports on natural resource quantities (such as forest area, cultivated land area, and change rate) and quality assessment reports (such as spatiotemporal changes in vegetation cover and soil moisture content); and third, trend results, including future natural resource change trend curves (predicted based on the ARIMA model) and risk warning reports (providing protection recommendations for areas with high change rates). The results are delivered in the form of GIS vector data, PDF reports, and visualization charts (line graphs, heat maps), meeting the survey, monitoring, and decision support needs of natural resource management departments.

[0103] The aforementioned integrated air-space-ground natural resource survey and monitoring method transforms digital elevation models into topological manifolds with navigational guidance value, turning terrain constraints, traditionally considered obstacles, into favorable conditions for UAV path planning. Based on satellite orbital mechanics prediction revisit windows, it establishes an integrated air-space-ground coupled optimization model. Distributed model predictive control enables collaborative task planning for UAV swarms. Within the limited time window of satellite overpasses, it constructs high-precision fingerprint maps using terrain features. Multi-UAV collaborative SLAM technology achieves relative positioning and stable formation maintenance in complex terrain environments where GNSS signals are completely denied. Finally, through deep fusion of multi-source air-space-ground data and temporal change detection and analysis, it effectively solves the technical bottlenecks of traditional survey methods in complex terrain environments, such as poor positioning accuracy, low data acquisition efficiency, and difficulties in air-space-ground coordination. This method represents a technological leap in natural resource survey and monitoring, from single-method to three-dimensional collaboration, from passive obstacle avoidance to active terrain utilization, and from decentralized operations to system optimization, significantly improving data acquisition quality and monitoring efficiency in complex environments.

[0104] refer to Figure 2 In one optional embodiment, a navigation topology manifold is constructed based on digital elevation model data of the target area, including the following steps:

[0105] S11. Based on digital elevation model data, ridgelines and valley lines are extracted through hydrological and topographic analysis and topographic curvature calculation, and then converted into topological graph nodes and topological graph edges.

[0106] Specifically, the digital elevation model (DEM) data is first preprocessed, and a priority flooding algorithm is used to fill the depressions in the DEM to avoid water flow "stagnation" in hydrological analysis. The filling threshold is set to the upper limit of the elevation error corresponding to the DEM resolution. At the same time, a bilateral filtering algorithm is used to remove random noise. This algorithm can preserve the terrain edge features while smoothing noise and prevent the outlines of ridge lines and valley lines from becoming blurred.

[0107] Valley lines were extracted through hydrotopographic analysis. The D8 algorithm was used to calculate the water flow direction for each pixel in the DEM, generating a water flow direction matrix. Based on this matrix, a recursive algorithm was used to calculate the cumulative runoff for each pixel, using the following formula: (in For pixels The cumulative amount of the confluence, (This is the set of all pixels upstream of the given pixel whose water flow direction points to that pixel). Set a threshold for the cumulative flow, and connect pixels with a cumulative flow greater than the threshold into a continuous line, which is the valley line. This corresponds to a low-lying channel and is suitable as a priority flight path for drones.

[0108] Ridge lines are extracted through topographic curvature calculation: Topographic curvature includes profile curvature and planar curvature. Profile curvature reflects the topographic convexity / concaveness along the water flow direction, while planar curvature reflects the topographic convexity / concaveness perpendicular to the water flow direction. The calculation formulas are as follows: , ,in , The slopes in the X and Y directions are calculated using the Sobel operator; , , The partial derivative is the second-order derivative, calculated using the finite difference method. The ridgeline identification condition is set as follows: profile curvature. (Convex slope) and plane curvature (Convex ridge), combined with slope threshold to exclude false identification of gentle areas, connect the pixels that meet the conditions into a continuous line, which is the ridge line, and its corresponding high terrain ridge line can be used as a reference boundary for UAV flight.

[0109] When converting ridgelines and valleylines into topological graph nodes and edges, the Douglas-Puk algorithm is used to simplify the lines, retaining inflection points, intersection points, and extreme points as core sampling points. The simplification threshold is set according to the DEM resolution to ensure that the geometric shape error of the simplified lines is small. These sampling points are defined as topological graph nodes, and each node is associated with its 3D coordinates in the DEM and its line type ("ridge node" or "valley node"). For any two nodes, if the planar distance is less than a preset neighborhood threshold (to ensure subsequent communication coverage) and there are no terrain obstacles in the connection, a topological graph edge is constructed, while duplicate edges and excessively short edges are deleted (to avoid frequent turning in path planning), forming the initial topological graph structure.

[0110] S12. Based on digital elevation model data, calculate the visibility index of the target area using a ray tracing algorithm; based on digital elevation model data, calculate the elevation change between different locations within the target area.

[0111] Specifically, the visibility index is calculated using a 3D ray tracing algorithm, for any pair of nodes in the topology graph. and Construct the ray parameter equations (plane projection: , Elevation direction: , Based on the DEM resolution, the ray sampling interval is determined to generate a set of sampling points. The actual terrain elevation of each sampling point is obtained from the DEM through bilinear interpolation. If all sampled points (excluding the node itself) satisfy... ( (A safety height threshold, set according to terrain complexity), then nodes are fully visible to each other; visibility index. If there is obstruction, Take the proportion of line-of-sight sampling points to the total number of sampling points. Range of values The higher the value, the better the visibility.

[0112] Elevation change calculation is performed on node pairs with already constructed edges in the topology graph, ensuring that the node elevation values ​​come from the same preprocessed DEM and have a consistent coordinate system. ( , (Node elevation). During the calculation, outlier elevations are checked. If the elevation exceeds the normal range for the target area, abnormal nodes are re-interpolated or removed. For node pairs with large planar distances, the effect of Earth's curvature is considered, and the formula is used to determine the appropriate value. ( For the Earth's radius, Correct the elevation value (for the horizontal distance between nodes), and then recalculate. .

[0113] S13. Based on the visibility index and combined with the elevation change, calculate the weight of each edge in the topological map to obtain the edge weight; the formula for calculating the edge weight is:

[0114]

[0115] in, This represents the edge weight of the topological graph edge corresponding to node i to j. The visibility index represents the path from node i to j in the topology graph. This represents the absolute value of the elevation change from node i to node j in the topology graph. This represents the horizontal projection distance from node i to node j in the topological graph. , , These are the weighting coefficients.

[0116] Specifically, For topology graph nodes arrive edge weights, For nodes arrive Visibility indicators For nodes arrive The absolute value of the elevation change, For nodes arrive The formula for the horizontal projection distance is: ).

[0117] S14. Construct a navigation topological manifold based on topological graph nodes, topological graph edges, and edge weights.

[0118] Specifically, the set of nodes and edges of the topological graph obtained in S11 are integrated with the edge weights calculated in S13 to form a weighted topological graph structure. During the integration process, feasible edges are selected based on their weights, and edges with excessively large weights (such as...) are removed. Paths exceeding preset thresholds (corresponding to extremely poor visibility, excessive undulation, or excessively long distances) retain edges whose weights align with the feasibility of UAV flight. Simultaneously, topological verification is performed on the retained edges and nodes to ensure that nodes form a connected network through feasible edges, avoiding isolated nodes (those unable to connect to other nodes via any feasible edges). The resulting navigation topology manifold transforms terrain constraints (such as visibility and elevation changes) into path cost information corresponding to edge weights, providing prior information for subsequent UAV swarm navigation and helping the swarm quickly select flight paths with strong terrain adaptability and low energy consumption.

[0119] In one alternative embodiment, a ternary coupled optimization model for the satellite-terrain-UAV swarm is established based on the revisit time window and ground coverage area, combined with the navigation topology manifold, including the following steps:

[0120] S21. With the objective of maximizing the coverage quality of the UAV swarm over the target area within the revisit time window, and considering the ground coverage range, an optimization objective function is constructed; the expression of the optimization objective function is:

[0121]

[0122] in, The total number of drones, Let be the start time of the revisit window for the k-th satellite. Let be the duration of the revisit time window for the k-th satellite. The spatial extent of the target area; For the i-th UAV at time t, the spatial point The coverage state function is constructed based on the ground coverage area; For spatial points The observation quality weighting function at time t.

[0123] Specifically, the logic of this optimization objective function is to ensure that the operations of the UAV swarm within the satellite revisit time window are deeply matched with the satellite's ground coverage, while prioritizing observation coverage of high-value areas. Among these, the coverage state function... The construction of associated satellite ground coverage: when space points When both are within the satellite's ground coverage area and within the sensor observation range of the i-th UAV (such as the camera's field of view or the lidar's scanning range), A value of 1 indicates that the point is effectively covered; if If the area is outside the satellite's ground coverage or outside the drone's observation range, then A value of 0 is used to prevent drones from operating ineffectively in areas without satellite coordination. Observation quality weighting function. The settings are based on the natural resource monitoring needs of the target area, such as canyon areas with complex terrain (high data collection difficulty and value) and key protected areas (high monitoring priority). Larger values ​​are taken for flat, bare land (low monitoring value). Smaller values ​​are used to prioritize coverage of key areas through weighted allocation.

[0124] S22. By adding constraints to the objective function, a ternary coupled optimization model is constructed; the constraints include UAV dynamics constraints, communication topology connectivity constraints, energy constraints, and satellite imaging matching constraints; the expression for the satellite imaging matching constraint is:

[0125]

[0126] in, For drones in location The observation angle at time t, For the satellite in position The observation angle at time t, This is the maximum permissible difference threshold between the observation angle of the UAV and the observation angle of the satellite.

[0127] Specifically, the purpose of adding constraints is to ensure that the optimization objectives meet the actual needs of UAV performance, communication coordination, energy reserves, and satellite imaging matching, forming a three-element coupled optimization model that balances "coverage quality" and "feasibility." UAV dynamics constraints are set based on the UAV's physical performance parameters, such as limiting flight speed to maximum speed, acceleration to maximum acceleration, and climb / descent angles to threshold values, preventing flight instability due to exceeding performance limits. Communication topology connectivity constraints, combined with navigation topology manifolds, require that the distance between any two adjacent UAVs does not exceed the maximum communication distance, ensuring the UAV swarm's communication network remains connected and supporting collaborative task data interaction. Energy constraints limit the total energy consumption of each UAV for the mission to no more than 80% of its remaining battery power (with 20% reserved for emergency return), and energy consumption calculations are based on flight path length and speed to prevent energy depletion mid-mission.

[0128] Satellite imaging matching constraints are the core constraints for ensuring the accuracy of air-space-ground data fusion. For drones in space The observation angle at time t usually refers to the angle between the optical axis of the UAV sensor and the direction of the ground normal (such as the camera tilt angle). The observation angle of the satellite at that point and time is calculated from the satellite's orbital parameters and the side-slip angle (e.g., the satellite's side-slip angle is...). At that time, ground point The observation angle is (Combined angle with orbital inclination) The maximum permissible difference threshold is set based on the data fusion accuracy requirements. This constraint limits the difference between the observation angles of the UAV and the satellite, preventing spatial misalignment between the acquired data and satellite imagery due to angular deviations, thus laying a geometric consistency foundation for subsequent air-space-ground data fusion.

[0129] In one optional embodiment, based on terrain fingerprinting, the UAV swarm is controlled to perform relative navigation and formation maintenance in a GNSS-denied environment through cooperative simultaneous localization and mapping (SMR), including the following steps:

[0130] S31. Based on the terrain fingerprint map, by establishing a backend optimization problem corresponding to the cooperative simultaneous localization and map building, a nonlinear least squares problem containing UAV pose constraints and map feature constraints is obtained.

[0131] Specifically, based on the constructed terrain fingerprint map, a corresponding backend optimization problem for collaborative simultaneous localization and map building is established. The terrain fingerprint map contains the spatial coordinates, local descriptors, and terrain attributes of terrain feature points. This information serves as the foundational data for backend optimization, establishing two core constraints: UAV pose constraints and map feature constraints. UAV pose constraints stem from the correlation of observations of the same terrain feature point by various UAVs at different times, ensuring that the pose changes of the same UAV within a continuous timeframe conform to motion patterns. Map feature constraints, on the other hand, are based on the consistency of observations of the same terrain feature point by multiple UAVs, ensuring the uniqueness of the terrain feature point's position in the global map. Integrating these two types of constraints transforms the problem into a nonlinear least squares problem involving UAV pose variables and map feature position variables. The core of this problem is to minimize the difference between the observed values ​​and the predicted values ​​calculated based on pose and map features.

[0132] S32. Solve the nonlinear least squares problem to obtain the optimized UAV pose estimate and map feature location estimate.

[0133] Specifically, the Gauss-Newton iterative algorithm or the Levenberg-Marquardt algorithm is used for solving the problem. During the iteration process, the estimated values ​​of the UAV pose variables and map feature position variables are continuously adjusted until the iteration error converges to a preset threshold (this threshold is set according to the navigation accuracy requirements, such as the position error converging to within 0.5 meters). Finally, the optimized UAV pose estimate (including the three-dimensional position and three-dimensional attitude angles of each UAV in the global coordinate system) and map feature position estimate (updating the accurate spatial coordinates of feature points in the terrain fingerprint map) are obtained, providing reliable basic data for subsequent navigation calculations.

[0134] S33. Based on UAV pose estimation and map feature position estimation, and combined with terrain constraints as navigation prior information, a height observation model is established by comparing the lidar ranging value with the digital elevation model data. Based on the height observation model, horizontal position constraints are established by using terrain feature lines to obtain the terrain-enhanced navigation solution.

[0135] Specifically, real-time ranging values ​​are acquired using the lidar onboard the UAV. These values ​​are then compared with the corresponding elevation values ​​in the digital elevation model (DEM) data for the target area. The deviation between these two values ​​is used to establish an altitude observation model, which corrects the UAV's vertical navigation errors. Simultaneously, horizontal position constraints are established using terrain feature lines (such as ridgelines and valley lines) in the terrain fingerprint map. The terrain feature lines observed by the UAV in real-time are matched with pre-stored feature lines in the map, and the matching deviation constrains the UAV's positional offset in the horizontal plane. Integrating the vertical constraints of the altitude observation model and the horizontal constraints of the terrain feature lines yields a terrain-enhanced navigation solution, which effectively compensates for the navigation accuracy degradation caused by missing GNSS signals.

[0136] S34. Perform spatiotemporal calibration on the navigation solution results to obtain a spatiotemporally consistent UAV swarm state estimate; based on the UAV swarm state estimate, generate UAV swarm formation maintenance control commands.

[0137] Specifically, in terms of time calibration, the unified time of the ground control station is used as the benchmark. Through the time synchronization protocol among the UAV swarm, the timestamp deviation of the navigation calculation results of each UAV is corrected to ensure that the state data of all UAVs are consistent in the time dimension. In terms of spatial calibration, based on the relative pose relationship between multiple UAVs, the spatial coordinate deviation caused by sensor errors of different UAVs is eliminated to obtain a spatiotemporally consistent UAV swarm state estimate (including the real-time position, velocity, and attitude of each UAV). Based on this state estimate, combined with preset formation geometry parameters (such as formation spacing and formation shape), the deviation between the expected pose and the actual pose of each UAV is calculated. Based on the deviation, formation maintenance control commands (such as velocity adjustment commands and attitude adjustment commands) are generated and sent to the flight control system of each UAV to achieve stable formation flight and relative navigation of the UAV swarm in GNSS denied environments.

[0138] In one optional embodiment, based on the terrain fingerprint map, a nonlinear least squares problem containing UAV pose constraints and map feature constraints is obtained by establishing a backend optimization problem corresponding to cooperative simultaneous localization and map building, including the following steps:

[0139] S41. Based on the location of map feature points in the terrain fingerprint map and the relative observation data between UAVs, define a state space including the UAV pose set and the map feature set.

[0140] Specifically, the location information of map feature points is first extracted from the terrain fingerprint map. This location information consists of the three-dimensional coordinates already marked in the map, which directly serve as the basic data for map features in the state space. Simultaneously, relative observation data is collected during the flight of the UAV swarm. This data is acquired by the lidar, millimeter-wave radar, or visual sensors carried by each UAV, covering key information such as the relative distance, relative angles (azimuth and pitch), and relative attitude angles between adjacent UAVs, reflecting the spatial relationship between them.

[0141] Based on the above two types of data, a state space is defined, which contains two core sets: one is the UAV pose set. ( (Number of drones), where the pose of each drone is... From three-dimensional position coordinates With three-dimensional attitude angle (roll angle) Pitch angle Yaw angle It consists of two parts: 1) a complete description of the UAV's position and attitude in space; and 2) a map feature set. ( (This refers to the number of map feature points), where each map feature point... Directly using the three-dimensional coordinates recorded in the topographic fingerprint map This ensures the consistency between map features and actual terrain features, providing an accurate reference benchmark for subsequent error calculations.

[0142] S42. Based on the UAV pose and map features in the state space, calculate the reprojection error of each UAV to each map feature point through the reprojection model, and calculate the relative pose error between UAVs through the relative pose model, so as to obtain the reprojection error set and the relative pose error set.

[0143] Specifically, for reprojection error calculation, a reprojection model is used to project the three-dimensional coordinates of map feature points onto the two-dimensional image plane of the UAV. Specifically, for the... The drone and the first Each map feature point, utilizing the drone pose The intrinsic parameters (focal length, principal point coordinates) and extrinsic parameters (pose parameters) are used to construct a reprojection model based on the principle of perspective projection. This model can calculate map feature points. In the The theoretical projected coordinates on the image plane of the drone. Compare these theoretical projected coordinates with the coordinates of the feature points actually observed by the drone. The reprojection errors are compared (i.e., the pixel coordinates of feature points extracted from the image), and the difference between the two is the single reprojection error. By traversing all drones and all map feature points, all calculated reprojection errors are summarized to form a reprojection error set.

[0144] For relative pose error calculation, a relative pose model is used to describe the pose relationship between adjacent UAVs. For the... frame and the first The drone, based on the pose of both and Constructing a relative pose model The model calculates the theoretical relative pose (including relative position and relative attitude) between two drones through coordinate transformation. This theoretical relative pose is then compared with the actual relative pose observed by the drones through sensors. The relative pose error is calculated by comparing the relative distance measured by lidar with the relative angle estimated by a visual sensor. The difference between the two is the single relative pose error. By iterating through all adjacent UAV pairs, all the calculated relative pose errors are summarized to form a set of relative pose errors.

[0145] S43. Based on the reprojection error set and the relative pose error set, and combining the covariance matrix corresponding to each error term, a nonlinear least squares problem is constructed. The expression for the nonlinear least squares problem is:

[0146]

[0147] in, , representing the set of pose states of all drones. Let N represent the pose state of the i-th drone, and N be the number of drones. , represents the set of all map feature points. This represents the three-dimensional coordinates of the k-th map feature point, where K is the number of map feature points; Let be the observation value of the i-th UAV at the k-th map feature point. Let i be the observation model function of the i-th UAV for the k-th map feature point. Let be the covariance matrix of the observations of the i-th UAV at the k-th map feature point. Let be the relative observation value between the i-th UAV and the (i+1)-th UAV. Let be the relative observation model function between the i-th UAV and the (i+1)-th UAV. Let be the covariance matrix of the relative observations between the i-th UAV and the (i+1)-th UAV; Represented by the covariance matrix The square of the defined Mahalanobis norm.

[0148] Specifically, This represents the set of pose states of all drones. For the first The pose state of the drone (including three-dimensional position and three-dimensional attitude). Represents the set of all map feature points. For the first The three-dimensional coordinates of a map feature point; It is the first The drone was used to attack the first The observed values ​​of map feature points (such as pixel coordinates in a two-dimensional image). It is the corresponding observation model function (i.e., the reprojection model). It is the covariance matrix of the observation, used to quantify the magnitude of observation noise (e.g., when visual observation noise is large). (The diagonal elements have larger values). It is the first frame and the first Relative observations between drones (such as relative distance and relative attitude). It is the corresponding relative observation model function. It is the covariance matrix of the relative observation, used to quantify the magnitude of the relative observation noise; Represented by the covariance matrix The square of the defined Mahalanobis norm is used to weight the error according to the magnitude of the observation noise. The larger the noise, the smaller the weight of the error term, thereby reducing the impact of noise on the optimization results.

[0149] The core objective of this nonlinear least squares problem is to adjust the drone pose set. With map feature set This minimizes the weighted sum of reprojection error and relative pose error, ultimately yielding an optimized result that balances UAV positioning accuracy and map feature accuracy, providing a reliable pose and map foundation for subsequent navigation calculations.

[0150] The aforementioned integrated air-space-ground natural resource survey and monitoring method transforms digital elevation models into topological manifolds with navigational guidance value, turning terrain constraints, traditionally considered obstacles, into favorable conditions for UAV path planning. Based on satellite orbital mechanics prediction revisit windows, it establishes an integrated air-space-ground coupled optimization model. Distributed model predictive control enables collaborative task planning for UAV swarms. Within the limited time window of satellite overpasses, it constructs high-precision fingerprint maps using terrain features. Multi-UAV collaborative SLAM technology achieves relative positioning and stable formation maintenance in complex terrain environments where GNSS signals are completely denied. Finally, through deep fusion of multi-source air-space-ground data and temporal change detection and analysis, it effectively solves the technical bottlenecks of traditional survey methods in complex terrain environments, such as poor positioning accuracy, low data acquisition efficiency, and difficulties in air-space-ground coordination. This method represents a technological leap in natural resource survey and monitoring, from single-method to three-dimensional collaboration, from passive obstacle avoidance to active terrain utilization, and from decentralized operations to system optimization, significantly improving data acquisition quality and monitoring efficiency in complex environments.

[0151] It should be understood that although the steps in the flowcharts of the embodiments described above 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 flowcharts of the embodiments described above may include multiple steps or multiple stages. These steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least some of the steps or stages of other steps.

[0152] Based on the same inventive concept, this application also provides an apparatus for implementing the aforementioned air-space-ground coordinated three-dimensional survey and monitoring method for natural resources. The solution provided by this apparatus is similar to the solution described in the above method; therefore, the specific limitations of one or more air-space-ground coordinated three-dimensional survey and monitoring apparatus embodiments provided below can be found in the above-described limitations of the air-space-ground coordinated three-dimensional survey and monitoring method for natural resources, and will not be repeated here.

[0153] In one exemplary embodiment, such as Figure 3 As shown, a three-dimensional natural resource survey and monitoring device 30 integrating air, space, and ground is provided to implement the methods in the above-described method embodiments. The device includes:

[0154] The terrain navigation manifold modeling module 31 is used to construct a navigation topology manifold based on the digital elevation model data of the target area; wherein, the navigation topology manifold is used to transform terrain constraints into prior information for UAV swarm navigation.

[0155] The satellite-terrain co-optimization module 32 is used to predict the satellite's revisit time window and ground coverage area based on the satellite orbital parameters and satellite imaging geometry characteristics using the Kepler orbital equations; and to establish a ternary coupled optimization model of satellite-terrain-UAV swarm based on the revisit time window and ground coverage area, combined with the navigation topology manifold.

[0156] The collaborative task planning module 33 solves the ternary coupled optimization model through a distributed model predictive control framework to generate a collaborative task planning scheme for the UAV swarm; the collaborative task planning scheme includes the waypoint sequence, velocity profile and task time window for each UAV.

[0157] The UAV formation self-reconstruction module 34 is used to control the UAV swarm to enter the target area based on the collaborative mission planning scheme within the revisit time window; and to perform topological self-reconstruction based on the navigation topology manifold to form an observation formation that is geometrically matched with satellite imaging.

[0158] The terrain feature extraction and navigation module 35 is used to extract terrain feature points in the target area and construct a terrain fingerprint map based on the terrain feature points. Based on the terrain fingerprint map, through coordinated simultaneous localization and map building, the module controls the UAV swarm to perform relative navigation and formation maintenance in a GNSS-denied environment.

[0159] The air-space-ground data fusion module 36 is used to control a swarm of drones to synchronously collect surface data based on the observation array when the satellite passes overhead, and to fuse the surface data with satellite imagery to obtain air-space-ground fused data; based on the change detection results and trend prediction results of the multi-period air-space-ground fused data, the results of the three-dimensional survey and monitoring of natural resources are generated.

[0160] Embodiments of this application also provide a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the aforementioned method embodiments.

[0161] Embodiments of this application also provide a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps in the above-described method embodiments.

[0162] For the device embodiments, since they basically correspond to the method embodiments, the relevant parts can be referred to in the description of the method embodiments. The device embodiments described above are merely illustrative. The components described as separate parts may or may not be physically separate, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this disclosure according to actual needs. Those skilled in the art can understand and implement this without creative effort.

[0163] The above-described embodiments are merely illustrative of several implementation methods of the embodiments of this application, and their descriptions are relatively specific and detailed. However, they should not be construed as limiting the scope of the patent application. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the embodiments of this application, and these modifications and improvements all fall within the protection scope of the embodiments of this application.

Claims

1. A three-dimensional survey and monitoring method for natural resources that integrates air, space, and ground, characterized in that, The method includes: S1. Construct a navigation topology manifold based on the digital elevation model data of the target area; wherein, the navigation topology manifold is used to transform terrain constraints into prior information for UAV swarm navigation; S2. Based on the satellite orbital parameters and satellite imaging geometry, predict the satellite's revisit time window and ground coverage area using the Kepler orbital equation; based on the revisit time window and the ground coverage area, establish a ternary coupled optimization model for the satellite-terrain-UAV swarm in conjunction with the navigation topology manifold; S3. Solve the ternary coupled optimization model using a distributed model predictive control framework to generate a collaborative task planning scheme for the UAV swarm; wherein, the collaborative task planning scheme includes the waypoint sequence, velocity profile, and task time window for each UAV; S4. Within the revisit time window, control the UAV swarm to enter the target area based on the cooperative mission planning scheme; perform topological self-reconstruction based on the navigation topology manifold to form an observation array that is geometrically matched with the satellite imaging; S5. Extract terrain feature points in the target area and construct a terrain fingerprint map based on the terrain feature points; based on the terrain fingerprint map, control the UAV swarm to perform relative navigation and formation maintenance in a GNSS denied environment through coordinated simultaneous localization and map building. S6. At the moment when the satellite passes overhead, the UAV swarm controlled by the observation array synchronously collects surface data, and fuses the surface data with the satellite imagery to obtain air-space-ground integrated data; based on the change detection results and trend prediction results of the air-space-ground integrated data over multiple periods, the results of the three-dimensional survey and monitoring of natural resources are generated.

2. The method according to claim 1, characterized in that, The navigation topology manifold is constructed based on the digital elevation model data of the target area, including: S11. Based on the digital elevation model data, extract ridge lines and valley lines through hydro-topographic analysis and topographic curvature calculation, and convert the ridge lines and valley lines into topological graph nodes and topological graph edges. S12. Based on the digital elevation model data, calculate the visibility index of the target area using a ray tracing algorithm; calculate the elevation change between different locations within the target area based on the digital elevation model data. S13. Based on the visibility index and the elevation change, calculate the weight of each edge in the topology map to obtain the edge weight; the formula for calculating the edge weight is: in, This represents the edge weight of the topological graph edge corresponding to node i to j. The visibility index represents the path from node i to j in the topology graph. This represents the absolute value of the elevation change from node i to node j in the topological graph. This represents the horizontal projection distance from node i to node j in the topology graph. , , These are the weighting coefficients; S14. Construct the navigation topology manifold based on the topology graph nodes, the topology graph edges, and the edge weights.

3. The method according to claim 1, characterized in that, The step of establishing a ternary coupled optimization model for satellite-terrain-UAV swarm based on the revisit time window and the ground coverage area, combined with the navigation topology manifold, includes: S21. With the objective of maximizing the coverage quality of the UAV swarm over the target area within the revisit time window, and in conjunction with the ground coverage range, construct an optimization objective function; the expression of the optimization objective function is: in, The total number of drones, Let be the start time of the revisit window for the k-th satellite. Let be the duration of the revisit time window for the k-th satellite. The spatial range of the target area; For the i-th UAV at time t, the spatial point A coverage state function, which is constructed based on the ground coverage area; For spatial points The observation quality weighting function at time t; S22. By adding constraints to the objective function, the ternary coupled optimization model is constructed; wherein, the constraints include UAV dynamics constraints, communication topology connectivity constraints, energy constraints, and satellite imaging matching constraints; the expression for the satellite imaging matching constraint is: in, For drones in location The observation angle at time t, For the satellite in position The observation angle at time t, This is the maximum permissible difference threshold between the observation angle of the UAV and the observation angle of the satellite.

4. The method according to any one of claims 1 to 3, characterized in that, Based on the terrain fingerprint map, the method of controlling a drone swarm to perform relative navigation and formation maintenance in a GNSS-denied environment through coordinated simultaneous localization and mapping includes: S31. Based on the terrain fingerprint map, by establishing a backend optimization problem corresponding to the cooperative simultaneous localization and the map construction, a nonlinear least squares problem containing UAV pose constraints and map feature constraints is obtained. S32. Solve the nonlinear least squares problem to obtain the optimized UAV pose estimate and map feature location estimate; S33. Based on the UAV pose estimation and the map feature position estimation, and combined with terrain constraints as navigation prior information, a height observation model is established by comparing the lidar ranging value with the digital elevation model data; based on the height observation model, a horizontal position constraint is established by using terrain feature lines to obtain the terrain-enhanced navigation solution. S34. Perform spatiotemporal calibration processing on the navigation solution results to obtain a spatiotemporally consistent UAV swarm state estimate; based on the UAV swarm state estimate, generate UAV swarm formation maintenance control commands.

5. The method according to claim 4, characterized in that, Based on the terrain fingerprint map, by establishing a backend optimization problem corresponding to the cooperative simultaneous localization and the map construction, a nonlinear least squares problem containing UAV pose constraints and map feature constraints is obtained, including: S41. Based on the map feature point locations in the terrain fingerprint map and the relative observation data between UAVs, define a state space including the UAV pose set and the map feature set; S42. Based on the UAV pose and map features in the state space, calculate the reprojection error of each UAV to each map feature point through the reprojection model, and calculate the relative pose error between UAVs through the relative pose model to obtain the reprojection error set and the relative pose error set. S43. Based on the reprojection error set and the relative pose error set, and combining the covariance matrix corresponding to each error term, construct the nonlinear least squares problem. The expression of the nonlinear least squares problem is: in, , representing the set of pose states of all drones. Let N represent the pose state of the i-th drone, and N be the number of drones. , represents the set of all map feature points. This represents the three-dimensional coordinates of the k-th map feature point, where K is the number of map feature points; Let be the observation value of the i-th UAV at the k-th map feature point. Let i be the observation model function of the i-th UAV for the k-th map feature point. Let be the covariance matrix of the observations of the i-th UAV at the k-th map feature point. Let be the relative observation value between the i-th UAV and the (i+1)-th UAV. Let be the relative observation model function between the i-th UAV and the (i+1)-th UAV. Let be the covariance matrix of the relative observations between the i-th UAV and the (i+1)-th UAV; Represented by the covariance matrix The square of the defined Mahalanobis norm.

6. A three-dimensional natural resource survey and monitoring device integrating air, space, and ground, for implementing the method described in any one of claims 1 to 5, characterized in that, The device includes: The terrain navigation manifold modeling module is used to construct a navigation topology manifold based on the digital elevation model data of the target area; wherein, the navigation topology manifold is used to transform terrain constraints into prior information for UAV swarm navigation; The satellite-terrain co-optimization module is used to predict the satellite's revisit time window and ground coverage area based on satellite orbital parameters and satellite imaging geometry characteristics using the Kepler orbital equations; and to establish a ternary coupled optimization model for the satellite-terrain-UAV swarm based on the revisit time window and the ground coverage area, combined with the navigation topology manifold. The collaborative task planning module solves the ternary coupled optimization model through a distributed model predictive control framework to generate a collaborative task planning scheme for the UAV swarm; wherein, the collaborative task planning scheme includes the waypoint sequence, velocity profile and task time window for each UAV; The UAV formation self-reconstruction module is used to control the UAV swarm to enter the target area based on the cooperative mission planning scheme within the revisit time window; and to perform topological self-reconstruction based on the navigation topology manifold to form an observation formation that is geometrically matched with the satellite imaging. The terrain feature extraction and navigation module is used to extract terrain feature points in the target area and construct a terrain fingerprint map based on the terrain feature points; based on the terrain fingerprint map, through coordinated simultaneous localization and map building, the module controls the UAV swarm to perform relative navigation and formation maintenance in a GNSS denied environment. The air-space-ground data fusion module is used to control a swarm of drones to synchronously collect surface data based on the observation array at the time of satellite overhead, and fuse the surface data with satellite imagery to obtain air-space-ground fused data; based on the change detection results and trend prediction results of the air-space-ground fused data over multiple periods, the module generates three-dimensional survey and monitoring results of natural resources.

7. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the method of any one of claims 1 to 5.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method of any one of claims 1 to 5.

Images