Three-dimensional terrain mapping system based on remote sensing technology

Abstract

The present invention provides a three-dimensional terrain mapping system based on remote sensing technology, which relates to the field of data processing technology, including: calculating the sun's zenith angle and azimuth according to the shooting time and geographical location information of remote sensing data; using the sun's zenith angle and azimuth to construct a second mapping relationship between a light source vector and the zenith angle and azimuth; calculating the reflected light intensity of each area according to the normal vector information and the light source vector in the first mapping relationship; obtaining the priori normal vector information of the surface of the target terrain area through inversion calculation based on the calculation of the reflected light intensity of each area and the remote sensing data; and correcting the normal vector information in the first mapping relationship based on the priori normal vector information to obtain the corrected normal vector information. The present invention can improve the accuracy and reliability of three-dimensional terrain mapping.

Description

Technical Field

[0001] The present invention relates to the technical field of data processing, in particular to a three-dimensional terrain mapping system based on remote sensing technology. Background Art

[0002] In passive remote sensing technology, reconstructing three-dimensional terrain by analyzing polarization information in remote sensing data is an effective method. Polarization information is generated when light waves interact with surface materials during propagation. It contains rich information about surface materials, such as surface roughness and dielectric constant, which is crucial for inferring surface morphology.

[0003] However, the existing three-dimensional terrain mapping system based on polarization information still has some defects in data processing:

[0004] For example, when extracting polarization data from remote sensing data, the extracted polarization data often contains a lot of noise and errors due to the influence of various factors such as atmospheric scattering and sensor noise. These noise and errors will directly affect the accuracy of subsequent terrain reconstruction.

[0005] For example, when constructing the mapping relationship between the normal vector and the incident light and light source vector, it is necessary to accurately estimate the incident angle and azimuth of the incident light, as well as the zenith angle and azimuth of the sun. However, due to the complexity of the shooting angle, time and geographic location information of remote sensing data, the estimation of these angles may have large errors, resulting in inaccurate mapping relationships. Summary of the invention

[0006] The technical problem to be solved by the present invention is to provide a three-dimensional terrain mapping system based on remote sensing technology, which can improve the accuracy and reliability of three-dimensional terrain mapping.

[0007] In order to solve the above technical problems, the technical solution of the present invention is as follows:

[0008] In a first aspect, a three-dimensional terrain mapping system based on remote sensing technology, the system comprises:

[0009] Data acquisition module, used to collect remote sensing data of different target terrain areas;

[0010] The data processing module is used to extract polarization data of the target terrain area from the remote sensing data; calculate the polarization degree of the terrain area surface by analyzing the polarization data; estimate the incident angle and azimuth of the incident light according to the shooting angle of the remote sensing data and the position of the sun; use the information of the polarization degree, the incident angle and the azimuth to construct a first mapping relationship between an initial normal vector and the incident angle and azimuth of the incident light; calculate the zenith angle and azimuth of the sun according to the shooting time and geographical location information of the remote sensing data; use the zenith angle and azimuth of the sun to construct a second mapping relationship between a light source vector and the zenith angle and azimuth; calculate the reflected light intensity of each area according to the normal vector information and the light source vector in the first mapping relationship; obtain the prior normal vector information of the target terrain area surface through inversion calculation based on the calculation of the reflected light intensity of each area and the remote sensing data; correct the normal vector information in the first mapping relationship based on the prior normal vector information to obtain the corrected normal vector information;

[0011] The three-dimensional reconstruction module is used to reconstruct the three-dimensional scene of the target area based on the corrected normal vector information and the polarization data of different target terrain areas at a preset angle.

[0012] Furthermore, remote sensing data of different target terrain areas are collected, including:

[0013] Divide the target terrain area into multiple sub-areas, each of which is responsible for one or more agents;

[0014] Assign initial remote sensing data collection parameters to each agent, including collection time, sensor configuration, and flight path;

[0015] Each agent simulates remote sensing data collection in the responsible sub-area according to its current collection parameters and evaluates the performance of the collection parameters to obtain local evaluation results;

[0016] Based on the information shared by all agents, a global estimation is performed to evaluate the overall performance of the current acquisition parameters to obtain a global estimation result;

[0017] According to the global estimation results, the acquisition parameters of each agent are updated, and the local estimation and global estimation are repeated until the preset number of iterations is reached to obtain an optimized remote sensing data acquisition plan;

[0018] According to the optimized remote sensing data collection plan, remote sensing data of different target terrain areas are collected.

[0019] Furthermore, the target terrain area is divided into multiple sub-areas, each of which is managed by one or more agents, including:

[0020] Determine the scope and boundaries of the target terrain area; preliminarily divide the target terrain area into multiple sub-areas based on the terrain's height, slope, vegetation cover and other characteristics, as well as the specific needs of remote sensing data collection;

[0021] For each initially divided sub-area, extract key terrain and ground features;

[0022] Analyze key terrain and landform features, and classify sub-areas according to the similarity of key terrain and landform features to obtain classification results;

[0023] According to the classification results, the spatial correlation between the various sub-regions and the continuity of the acquisition tasks are analyzed to obtain the sub-region groups with correlation;

[0024] Each sub-region that has been classified and determined to be relevant is regarded as a node in the graph model, and the feasible paths between sub-regions are regarded as edges in the graph;

[0025] Define the rules for ants to move in the graph, initialize the ant population, and each ant chooses the next sub-area to move to based on the pheromone concentration and heuristic information of the nodes around the current location. Ants release pheromones when selecting a path and update the pheromone concentration based on the length and quality of the path.

[0026] When all ants have completed a path search, they will globally update the pheromone according to the length and quality of the overall path until the preset number of iterations is reached. After the iteration is completed, the optimized sub-area division scheme and the corresponding path planning results will be output;

[0027] According to the optimized sub-area division scheme and the corresponding path planning results, the intelligent agent is assigned to the corresponding sub-area for remote sensing data collection.

[0028] Furthermore, the key terrain and landform features are analyzed, and the sub-regions are classified according to the similarity of the key terrain and landform features to obtain classification results, including:

[0029] The terrain and object features of each sub-region are combined into a feature vector, which represents the comprehensive features of the sub-region;

[0030] Determine the corresponding K value through cross-validation, prepare an empty result list for the sub-regions to be classified; for each sub-region to be classified, perform the following steps:

[0031] Calculate the distance between the feature vector of the current sub-region and the feature vectors of all sub-regions in the training set;

[0032] Sort the distances and select the corresponding K sub-regions as the neighbors of the current sub-region to obtain the classification result.

[0033] Furthermore, according to the classification results, the spatial correlation between the various sub-regions and the continuity of the acquisition tasks are analyzed to obtain sub-region groups with correlation, including:

[0034] Analyze the adjacency relationship between each sub-region, and for sub-regions that are not directly adjacent, calculate the spatial distance between the sub-regions;

[0035] Based on the adjacency relationship and spatial distance, a topological relationship network between sub-regions is established, and the topological relationship network represents the spatial connectivity of sub-regions;

[0036] According to the goal of the collection task, evaluate the continuity of task execution between different sub-areas to obtain the analysis results of the continuity of the collection task;

[0037] Based on the analysis results of the topological relationship network and the continuity of the acquisition task, the sub-areas are comprehensively evaluated to obtain a comprehensive evaluation result;

[0038] Based on the comprehensive evaluation results, a grouping strategy is formulated to obtain relevant sub-region groups.

[0039] Furthermore, polarization data of the target terrain area is extracted from the remote sensing data; and the polarization degree of the terrain area surface is calculated by analyzing the polarization data, including:

[0040] Obtaining raw data containing polarization information from remote sensing satellites or ground-based remote sensing equipment;

[0041] Filtering the dataset covering the target terrain area from the original data;

[0042] According to the characteristics of remote sensing data, the polarization band containing polarization information is determined;

[0043] For each polarization band, the corresponding polarization component is extracted, and the polarization bands are combined into a polarization data set;

[0044] For each region, calculate its degree of polarization;

[0045] The polarization degree is assigned to the corresponding area to generate a polarization degree map.

[0046] Further, according to the normal vector information and the light source vector in the first mapping relationship, the reflected light intensity of each area is calculated, including:

[0047] Obtain normal vector information of each region from the first mapping relationship, and calculate the light source vector according to the zenith angle and azimuth angle of the sun;

[0048] According to the normal vector information in the first mapping relationship and the calculated light source vector, Calculate the reflected light intensity for each area ;

[0049] in, Indicates the direction of incident light; Indicates the viewing direction; represents the normal vector; represents the roughness parameter; is the half vector of incident and reflected light; ; is a parameter related to surface roughness and its value is ;in, ;in, is the reflectivity at normal incidence.

[0050] Furthermore, based on the calculation of the reflected light intensity and remote sensing data of each area, the prior normal vector information of the surface of the target terrain area is obtained through inversion calculation, including:

[0051] Obtain the surface reflectance data of the target terrain area, and calculate the relationship between the incident angle and the reflection angle based on the sun position, sensor viewing angle and surface geometry;

[0052] The residual is defined as the difference between the observed reflected light intensity and the light intensity predicted by the parameterized reflectance process;

[0053] Construct the objective function of the residual sum of squares and set the initial values ​​for the unknown parameters in the reflection process;

[0054] The gradient descent method is used to minimize the residual sum of squares through iterative calculation. In each iteration, the parameter value is adjusted;

[0055] The residual sum of squares between the predicted data and the observed data after each iteration is calculated. When the change in the residual sum of squares reaches the preset number of iterations, the iteration is stopped and the parameter value in the current iteration step is output as the final solution, where the final solution of the normal vector is the prior normal vector information of the surface of the target terrain area.

[0056] Furthermore, the objective function The calculation formula is:

[0057] ;

[0058] in, is the number of observation points; At the observation point The observed intensity of reflected light; is the weight of each observation point; It is an observation point The standard deviation of the observation error; is the regularization parameter; is the angle between the incident light and the surface normal; is the azimuth angle between the incident light and the reflected light; Represents the index value; is the reflectivity of the surface.

[0059] In the second aspect, a three-dimensional terrain mapping method based on remote sensing technology is applied to the system, including:

[0060] Collect remote sensing data of different target terrain areas;

[0061] Extract polarization data of the target terrain area from remote sensing data; calculate the polarization degree of the terrain area surface by analyzing the polarization data; estimate the incident angle and azimuth of the incident light according to the shooting angle of the remote sensing data and the position of the sun; use the information of the polarization degree, incident angle and azimuth to construct a first mapping relationship between an initial normal vector and the incident angle and azimuth of the incident light; calculate the zenith angle and azimuth of the sun according to the shooting time and geographical location information of the remote sensing data; use the zenith angle and azimuth of the sun to construct a second mapping relationship between a light source vector and the zenith angle and azimuth; calculate the reflected light intensity of each area according to the normal vector information and the light source vector in the first mapping relationship; obtain the prior normal vector information of the surface of the target terrain area through inversion calculation based on the calculation of the reflected light intensity of each area and the remote sensing data; correct the normal vector information in the first mapping relationship based on the prior normal vector information to obtain the corrected normal vector information;

[0062] The three-dimensional scene of the target area is reconstructed based on the corrected normal vector information and the polarization data of different target terrain areas at preset angles.

[0063] The above scheme of the present invention includes at least the following beneficial effects.

[0064] Through the data acquisition module, the system can quickly and efficiently collect remote sensing data of different target terrain areas. This automated data acquisition method not only reduces the complexity and cost of manual operations, but also greatly improves the speed and scope of data acquisition.

[0065] The data processing module extracts polarization data from remote sensing data through refined steps and calculates the polarization degree of the terrain area. Furthermore, by constructing the mapping relationship between the normal vector and the incident light and light source vector, and performing inversion calculations using the reflected light intensity and remote sensing data, the system can obtain more accurate prior normal vector information and correct the normal vector. This series of processing procedures significantly improves the accuracy and reliability of data processing.

[0066] The 3D reconstruction module uses the corrected normal vector information and polarization data at a preset angle to reconstruct the 3D scene of the target area. Due to the accuracy of the previous data processing, the 3D reconstruction module can generate a more realistic and detailed terrain model that accurately reflects the undulations, textures and other detailed features of the target terrain. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 It is a schematic diagram of a three-dimensional terrain mapping system based on remote sensing technology provided by an embodiment of the present invention.

[0068] Figure 2 It is a schematic flow chart of a three-dimensional terrain surveying and mapping method based on remote sensing technology provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0069] The exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present disclosure are shown in the accompanying drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided to enable a more thorough understanding of the present disclosure and to fully convey the scope of the present disclosure to those skilled in the art.

[0070] like Figure 1 As shown, the embodiment of the present invention proposes a three-dimensional terrain mapping system based on remote sensing technology, including:

[0071] Data acquisition module, used to collect remote sensing data of different target terrain areas;

[0072] The data processing module is used to extract polarization data of the target terrain area from the remote sensing data; calculate the polarization degree of the terrain area surface by analyzing the polarization data; estimate the incident angle and azimuth of the incident light according to the shooting angle of the remote sensing data and the position of the sun; use the information of the polarization degree, the incident angle and the azimuth to construct a first mapping relationship between an initial normal vector and the incident angle and azimuth of the incident light; calculate the zenith angle and azimuth of the sun according to the shooting time and geographical location information of the remote sensing data; use the zenith angle and azimuth of the sun to construct a second mapping relationship between a light source vector and the zenith angle and azimuth; calculate the reflected light intensity of each area according to the normal vector information and the light source vector in the first mapping relationship; obtain the prior normal vector information of the target terrain area surface through inversion calculation based on the calculation of the reflected light intensity of each area and the remote sensing data; correct the normal vector information in the first mapping relationship based on the prior normal vector information to obtain the corrected normal vector information;

[0073] The three-dimensional reconstruction module is used to reconstruct the three-dimensional scene of the target area based on the corrected normal vector information and the polarization data of different target terrain areas at a preset angle.

[0074] In an embodiment of the present invention, the data acquisition module can quickly capture and collect remote sensing data of the target terrain area, effectively reducing the time cost of data acquisition and improving the efficiency of the overall mapping work. The module can adapt to the remote sensing data acquisition needs of different terrain areas, whether it is a city, mountainous area or water area, it can effectively capture data, showing extremely high adaptability and flexibility. By finely processing the polarization information in the remote sensing data, the data processing module can accurately calculate the polarization degree of the terrain surface, providing key parameters for constructing an accurate normal vector mapping relationship, thereby significantly improving the accuracy of three-dimensional terrain mapping. The module successfully constructs the mapping relationship between the normal vector and the incident light and light source vector. This innovative method enhances the scientificity and accuracy of data processing. Based on the prior normal vector information obtained by inversion calculation, the initial normal vector is corrected, which further improves the accuracy of the normal vector and ensures the reliability and authenticity of the three-dimensional reconstruction. Using the corrected normal vector information and polarization data, the three-dimensional reconstruction module can generate a highly realistic three-dimensional terrain model, carefully restore the undulations, textures and lighting effects of the terrain, and greatly enhance the realism of the visual experience.

[0075] In a preferred embodiment of the present invention, collecting remote sensing data of different target terrain areas includes:

[0076] Divide the target terrain area into multiple sub-areas, each of which is responsible for one or more agents;

[0077] Assign initial remote sensing data collection parameters to each agent, including collection time, sensor configuration, and flight path, including: determine the scope and boundaries of each sub-area, and then assign one or more agents to each sub-area to ensure that each sub-area has a responsible agent; set an initial collection time window for each agent based on the geographic location, season, and weather conditions of the target terrain area to ensure data collection under optimal lighting and meteorological conditions; select appropriate sensor configurations for each agent based on the required data type (such as multispectral, hyperspectral, radar, etc.) and resolution requirements, including sensor type, band settings, and calibration parameters; plan an initial flight path for each agent based on the terrain characteristics of the sub-area, obstacle distribution, and sensor coverage to ensure path safety and data collection integrity;

[0078] Each agent simulates remote sensing data acquisition in the responsible sub-area according to its current acquisition parameters, and evaluates the performance of the acquisition parameters to obtain local evaluation results, which specifically include: each agent simulates remote sensing data acquisition in the responsible sub-area according to its current acquisition parameters (acquisition time, sensor configuration and flight path) to generate simulated data; the quality of the simulated data is evaluated, including indicators such as data integrity, clarity and accuracy, to evaluate the performance of the current acquisition parameters in the local range; the local evaluation results of each agent are summarized and fed back to the global estimation module; the local evaluation of the agent can be specifically calculated by the following formula:

[0079] ;

[0080] in, Indicates the effective data collection area; Indicates the effective scheduled data collection area; Indicates The maximum brightness value in a pixel neighborhood; Indicates The minimum brightness value in a pixel neighborhood; Represents the total number of pixels in the image; represents the root mean square error between the simulated data and the real observed data; represents the Pearson correlation coefficient between the simulated data and the observed data; , , and represents the weight coefficient; Represents an index; Represents the local evaluation score of the agent.

[0081] Based on the information shared by all agents, a global estimation is performed to evaluate the overall performance of the current acquisition parameters to obtain a global estimation result, which specifically includes: collecting the local evaluation results shared by all agents, as well as other possible relevant information (such as terrain complexity, weather forecast, etc.); based on the integrated information, a global performance evaluation is performed, taking into account the efficiency, quality and cost of data acquisition, and obtaining the overall performance evaluation result of the current acquisition parameters; the global estimation result is output to each agent to provide a basis for the next step of parameter update; the global estimation result can be specifically calculated by the following calculation formula:

[0082] ;

[0083] in, is the number of agents; Indicates The local evaluation score of each agent; Represents the global estimation result.

[0084] According to the global estimation results, the acquisition parameters of each agent are updated, and the local estimation and global estimation are repeated until the preset number of iterations is reached to obtain an optimized remote sensing data acquisition plan, which includes: by analyzing the global estimation results, the key issues affecting the overall performance are identified, such as poor data acquisition quality in some areas, inefficient flight paths, or unreasonable sensor configuration; according to the problems identified in the global estimation results, the adjustment goals are determined, such as improving data acquisition quality, optimizing flight efficiency, or reducing costs; for each adjustment goal, the acquisition parameters to be adjusted are selected, and the parameters that may be involved include acquisition time window, sensor configuration (such as resolution, band selection, etc.) and flight path. According to the feedback of the global estimation results, the adjustment direction and amplitude of each parameter are determined. For example, if the acquisition quality of a certain area is found to be poor, the sensor configuration needs to be adjusted to improve the resolution, or the flight path needs to be optimized to ensure better coverage. According to the global estimation results, the acquisition time is adjusted to avoid adverse meteorological conditions or improve the lighting conditions for data acquisition. According to the identified problems, the configuration parameters of the sensor are adjusted, such as changing the resolution, adjusting the band settings, or optimizing the calibration parameters, to improve the quality and accuracy of data acquisition. Based on the feedback from the global estimate, the flight path is optimized to improve the collection efficiency and coverage. This may involve replanning the flight path, adjusting the flight altitude and speed, etc. The adjusted collection parameters are applied to each agent to ensure that each agent collects data according to the new parameters.

[0085] After updating the parameters, the iterative process of local estimation and global estimation is repeated. Each agent collects simulated data in the sub-region it is responsible for and evaluates the performance of the new parameters. Then, a global estimation is performed based on the shared information of all agents. The iterative process is repeated until the preset number of iterations is reached or a specific stopping condition is met, such as the global performance indicator reaches a preset threshold or the performance improvement of several consecutive iterations is less than a certain threshold.

[0086] According to the optimized remote sensing data collection plan, remote sensing data of different target terrain areas are collected, including: after multiple iterations of optimization, the final remote sensing data collection plan is obtained, including detailed information such as the specific collection time, sensor configuration and flight path of each intelligent agent; according to the final collection plan, each intelligent agent performs actual remote sensing data collection work according to the specified time and path to ensure the accuracy and completeness of the data.

[0087] In an embodiment of the present invention, by dividing the target terrain area into multiple sub-areas and assigning an agent to each sub-area for data collection, the method realizes parallel processing and greatly improves the efficiency of data collection. At the same time, by iteratively optimizing the acquisition parameters, it can ensure that the final remote sensing data acquisition plan is carefully optimized, so as to achieve the best balance between resource utilization and data quality. The method assigns initial remote sensing data acquisition parameters to each agent. This flexibility enables the system to make adaptive adjustments according to the characteristics and real-time conditions of different terrain areas, thereby ensuring the pertinence and effectiveness of data collection. Through the mechanism of combining local evaluation and global estimation, the method can comprehensively evaluate the performance of the current acquisition parameters, timely discover and correct possible problems, and this mechanism helps to improve the quality of the remote sensing data finally collected; the optimized remote sensing data acquisition plan can more reasonably allocate resources, such as sensor configuration and flight path, so as to effectively reduce resource consumption and cost expenditure while ensuring data quality.

[0088] In a preferred embodiment of the present invention, the target terrain area is divided into multiple sub-areas, each of which is managed by one or more agents, including:

[0089] Determine the scope and boundaries of the target terrain area; preliminarily divide the target terrain area into multiple sub-areas based on the terrain's height, slope, vegetation cover and other characteristics, as well as the specific needs of remote sensing data collection, including: Collect relevant data on the target terrain area, which usually includes map data, satellite images or geographic information system (GIS) data; Use the collected data to identify the specific scope and boundaries of the target terrain area in GIS software or special mapping tools, which can be done by drawing polygons, lines or using existing geographic boundary data; Assign geographic coordinates (such as longitude and latitude) to each key point of the boundary to ensure the accuracy of the boundary; Use remote sensing technology (such as satellite remote sensing, drone remote sensing, etc.) to obtain detailed terrain data such as height, slope, vegetation cover, etc. of the target terrain area, which are represented in the form of digital elevation model (DEM), slope map and vegetation index map; Determine the specific needs of remote sensing data collection, including the type, resolution, timeliness, etc. of the required data; Combined with the terrain characteristics and collection requirements, use the spatial analysis tools in GIS software (such as cluster analysis, overlay analysis, etc.) to preliminarily divide the target terrain area. For example, the area can be divided into different elevation zones according to the height range, or steep areas can be separated from flat areas according to the slope threshold; based on the preliminary division, the boundaries of the sub-areas are further refined to ensure that each sub-area has relative consistency in terrain characteristics, which can be achieved by smoothing boundaries, merging similar areas, or splitting complex areas.

[0090] For each initially divided sub-area, extract key terrain and landform features, including: determine the key terrain and landform features to be extracted, including average elevation, maximum or minimum slope, main vegetation types and water body distribution, among which the average elevation is used to understand the overall elevation of the sub-area; the maximum or minimum slope reflects the steepness of the terrain, affecting the difficulty and safety of data collection; the main vegetation types are used to identify forests, grasslands, farmlands, etc., which are crucial for surface cover and ecological environment analysis; water body distribution includes rivers, lakes, reservoirs, etc.; according to the specific purpose of remote sensing data collection (such as environmental monitoring, resource survey, urban planning, etc.); convert the pixel values ​​of the remote sensing image into physically meaningful radiation brightness values ​​or reflectance values, using the calibration coefficients or lookup tables provided by the sensor for conversion; eliminate the influence of atmospheric scattering and absorption on the remote sensing image to obtain the true reflectance of the surface; use the atmospheric radiation transmission model (such as M ODTRAN, 6S, etc.) or image-based atmospheric correction methods (such as dark pixel method, histogram matching method, etc.); correct the geometric distortion of remote sensing images caused by factors such as sensor attitude, earth rotation, and terrain undulation; use ground control points (GCPs) and geometric transformation models (such as polynomial model, rational function model, etc.) for correction; select training samples and use classifiers (such as maximum likelihood method, support vector machine, etc.) to classify remote sensing images to obtain vegetation type maps; use digital elevation model (DEM) data to calculate the average elevation through the terrain analysis tool of GIS software; use the same DEM data to calculate the maximum or minimum slope through the slope analysis tool; use the spectral characteristics of water bodies in remote sensing images (such as low reflectivity, high absorption rate, etc.) to extract water body boundaries through threshold segmentation or shape analysis, and combine DEM data to further identify and analyze the three-dimensional morphology and distribution characteristics of water bodies.

[0091] Analyze key terrain and landform features, and classify sub-areas according to the similarity of key terrain and landform features to obtain classification results;

[0092] According to the classification results, the spatial correlation between the various sub-regions and the continuity of the acquisition tasks are analyzed to obtain the sub-region groups with correlation;

[0093] Each sub-region determined by classification and correlation is regarded as a node in the graph model, and the feasible paths between sub-regions are regarded as edges in the graph, specifically including: each sub-region determined by classification and correlation is regarded as a node in the graph, and each node stores its corresponding sub-region information, such as location, characteristics, etc.; if there is a feasible path between two sub-regions (for example, geographically adjacent and without obstacles), an edge is added between the corresponding nodes, and the edge can have a weight to indicate the length of the path; according to the definition of nodes and edges, a complete graph model is constructed to represent the topological structure of the entire region.

[0094] Define the rules for ants to move in the graph, initialize the ant population, and each ant selects the next sub-area to move to based on the pheromone concentration and heuristic information of the nodes around the current position. The ants release pheromones when selecting a path and update the pheromone concentration based on the length and quality of the path. Specifically, it includes: initializing a certain concentration of pheromones on each edge in the graph. The pheromone concentration can represent the attractiveness of the path; calculating heuristic information for each node, such as the straight-line distance to the target node or the estimated movement cost based on the terrain features, randomly selecting the starting node in the graph, and initializing a certain number of ants. Each ant has a memory ability to record the path it has traveled. Each ant uses a probabilistic selection formula to select the next node to move based on the pheromone concentration and heuristic information of the nodes around the current position. Among them, after the ant path search, the update of the pheromone concentration on each edge can be expressed by the following formula:

[0095] ;

[0096] in, Indicates that at time step When, the path pheromone concentration on; Indicates pheromone volatility ;in, Indicates the path The amount of pheromone released can be calculated based on the length and quality of the path:

[0097] ;

[0098] in, is a constant, representing the total amount of pheromone intensity; represents the path length or movement cost, where the calculation formula for pheromone intensity is: , represents the coefficient, Indicates the path pheromone concentration on; Indicates the path The pheromone strength of . The probability of each ant choosing the next node can be expressed by the following formula:

[0099] ;

[0100] in, Indicates ant slave node Move to Node The probability of Indicates the path pheromone concentration on; Indicates the path pheromone concentration on; Represents heuristic information; Represents heuristic information, which can be represented by the distance between nodes or the movement cost, and is usually defined as:

[0101] ;

[0102] in, Representation Node and nodes The distance between Representation Node The set of neighbor nodes of and Represents the coefficient.

[0103] When all ants complete a path search, the pheromone is globally updated according to the length and quality of the overall path until the preset number of iterations is reached. After the iteration, the optimized sub-area division scheme and the corresponding path planning results are output, which specifically include: each ant searches for a path in the graph according to the defined movement rules until it reaches the target node or meets other termination conditions. During the movement of the ants, pheromones are released on the edges they pass through according to the length and quality of the path (such as movement cost, time, etc.), and the pheromone concentration is locally updated. Shorter paths or paths with higher quality should obtain more pheromones. When all ants complete a path search, the pheromone is globally updated for all edges in the graph according to the length and quality of the overall path, which involves the volatilization of pheromones and enhancement according to the found optimal path. Repeat the steps of ant path search and pheromone update until the preset number of iterations is reached or other stopping criteria are met (such as finding the optimal solution that meets the requirements). After the iteration, the sub-area division scheme is determined according to the optimal path found by the ants, and the corresponding path planning results are output, which may include information such as the boundary, center point, and connection path of each sub-area.

[0104] According to the optimized sub-area division scheme and the corresponding path planning results, the intelligent agents are assigned to the corresponding sub-areas for remote sensing data collection, specifically including: according to the optimized sub-area division scheme and the corresponding path planning results, the intelligent agents (such as drones, ground robots, etc.) are assigned to the corresponding sub-areas for remote sensing data collection, so that each sub-area is covered by a suitable intelligent agent. The intelligent agents collect remote sensing data in their respective sub-areas according to the planned path, and collect the required terrain, ground objects and other information. During the data collection process, the status and collection progress of the intelligent agents are monitored in real time, and necessary adjustments are made according to the actual situation, such as re-planning the path, reinforcing the intelligent agents, etc., to ensure the smooth completion of the collection task.

[0105] In the embodiment of the present invention, through reasonable sub-area division, it can ensure that the area responsible for each agent is more concentrated and continuous, reduce the invalid movement of agents between different sub-areas, and thus improve the overall data collection efficiency. By classifying the sub-areas according to the similarity of terrain and object features, agents can be more targeted and collection parameters can be adjusted, so as to achieve optimal allocation of resources and avoid waste of resources. By analyzing the key terrain and object features, areas that need to be focused on can be more accurately identified, so that more refined data collection can be performed in these areas, thereby improving the overall data quality. The application of the ant colony algorithm makes path planning more intelligent and flexible, and can adjust the path selection according to real-time conditions, avoiding the limitations that may be brought about by fixed path planning, and improving the rationality and adaptability of path planning. The optimized path planning reduces the moving distance and time of the agent, thereby reducing energy consumption and maintenance costs; by dividing the target area into multiple sub-areas and assigning an independent agent to each sub-area, the redundancy and fault tolerance of the system are increased. Even if an agent fails, it will not affect the operation of the entire system.

[0106] In a preferred embodiment of the present invention, key terrain and landform features are analyzed, and sub-regions are classified according to the similarity of key terrain and landform features to obtain classification results, including:

[0107] The terrain and landform features of each sub-region are combined into a feature vector to represent the comprehensive features of the sub-region, including: for each sub-region, collect relevant terrain and landform feature data, which may include elevation, slope, surface cover type (such as forest, water body, bare land, etc.), building density, etc.; extract key indicators that can represent the characteristics of the sub-region from the collected data, for example, the average elevation and maximum slope in the sub-region can be calculated by statistical methods; combine the extracted features into a feature vector. This vector can be multi-dimensional, with each dimension representing a feature. For example, a feature vector may contain elements such as [average elevation, maximum slope, forest coverage, building density].

[0108] Determine the corresponding K value through cross-validation, prepare an empty result list for the sub-region to be classified; for each sub-region to be classified, perform the following steps: calculate the distance between the feature vector of the current sub-region and the feature vectors of all sub-regions in the training set; sort the distances, select the corresponding K sub-regions as the neighbors corresponding to the current sub-region, to obtain the classification result, specifically including:

[0109] Divide the sub-region data set of known classification into training set and validation set, which can be done by random division or stratified sampling. According to experience, set a search range of K value, such as K = 1 to K = 10. For each K value, use the training set to train the K nearest neighbor classifier, and evaluate the performance of the classifier on the validation set. The performance evaluation indicators can be accuracy, recall rate, F1 score, etc. According to the results of cross-validation, select the K value with the best performance as the parameter of the K nearest neighbor algorithm in the subsequent classification task. Create an empty result list for the sub-region to be classified to store the classification results of each sub-region.

[0110] For each sub-region to be classified, perform the following operations:

[0111] Get the feature vector of the current sub-region, and use the Euclidean distance to calculate the distance between the feature vector of the current sub-region and the feature vectors of all sub-regions in the training set; sort the calculated distances from small to large to obtain the sorting result; according to the sorting result, select the K sub-regions with the smallest distance as the neighbors of the current sub-region to be classified; according to the known classification of these K neighbors, determine the classification result of the current sub-region to be classified by voting or other decision rules (such as weighted voting); add the classification result to the previously prepared empty result list; when all sub-regions to be classified have been classified, sort the result list to ensure that the classification result of each sub-region corresponds to it.

[0112] In an embodiment of the present invention, through a detailed analysis of the terrain and land features, the similarities and differences between different sub-regions can be captured more accurately, thereby achieving more accurate classification. This classification method reduces the possibility of misclassification and improves the accuracy of the overall classification. As a representative of the comprehensive characteristics of the sub-region, the feature vector helps to explain the classification results, making the classification process more transparent and credible; accurate classification helps decision makers better understand the characteristics and needs of different sub-regions, thereby optimizing resource allocation. For example, in the fields of environmental protection, urban planning, disaster prevention, etc., more targeted strategies and measures can be formulated based on the terrain and land features of different sub-regions. By automatically processing large amounts of data, the method can quickly provide classification results for each sub-region, thereby speeding up the decision-making process.

[0113] In a preferred embodiment of the present invention, according to the classification results, the spatial correlation between the various sub-regions and the continuity of the acquisition tasks are analyzed to obtain a sub-region group with correlation, including:

[0114] Analyze the adjacency relationship between each sub-region, and for non-directly adjacent sub-regions, calculate the spatial distance between the sub-regions, specifically including: obtaining the geographical location and boundary information of all sub-regions and storing them in GIS (Geographic Information System) data or spatial database; using GIS tools or spatial analysis algorithms to detect whether each sub-region is directly adjacent to other sub-regions, which can be achieved by checking whether the sub-region boundaries share edges or points; for each pair of directly adjacent sub-regions, record their adjacency relationship, which can be an adjacency matrix or an adjacency list; for each pair of non-directly adjacent sub-regions, use Euclidean distance to calculate the shortest spatial distance between them; and store the calculated distance data in a distance matrix.

[0115] Based on the adjacency relationship and spatial distance, a topological relationship network is established between sub-regions. The topological relationship network represents the spatial connectivity of sub-regions, which includes: creating an empty graph or network data structure to represent the topological relationship between sub-regions; adding each sub-region as a node in the network to the data structure; adding edges to the network based on the adjacency relationship and the calculated spatial distance. Undirected edges are added between directly adjacent sub-regions, and the weight can be set to 1 or 0; for non-directly adjacent sub-regions, directed edges or undirected edges can be added based on the distance, and the weight can be set to the distance value between them; check whether the constructed network correctly represents the spatial connectivity between sub-regions.

[0116] According to the objectives of the acquisition task, the continuity of task execution between different sub-areas is evaluated to obtain the analysis results of the continuity of the acquisition task, including: evaluating the similarity of the two sub-areas in terrain, such as elevation changes, slopes, etc., terrain data (such as DEM digital elevation model) can be used to calculate these indicators; comparing whether the types of land objects in the two sub-areas are consistent, such as water bodies, vegetation, buildings, etc., which can be obtained through remote sensing image classification or land use data; evaluating the shortest path or travel time from one sub-area to another, for example, GIS network analysis tools can be used to calculate the shortest path or travel time; collecting and organizing the required GIS data, including terrain data, remote sensing images, land use data, transportation network data, etc.; selecting appropriate GIS analysis methods according to the defined indicators. For example, use spatial interpolation or terrain analysis tools to calculate terrain similarity; use remote sensing image classification technology to extract land features; for each pair of adjacent or similarly distant sub-regions, calculate their continuity index values ​​separately; organize the calculated continuity index values ​​into a standard format; based on the calculation results, construct a continuity matrix or list, in which each element represents the continuity value between a pair of sub-regions, which can be in the form of a two-dimensional array or table for easy storage and query.

[0117] Based on the analysis results of the topological relationship network and the continuity of the acquisition task, a comprehensive evaluation is conducted on the sub-areas to obtain a comprehensive evaluation result, including: making the data format of the topological relationship network data and the task continuity analysis results consistent; associating each node (sub-area) in the topological relationship network with the corresponding sub-area in the task continuity analysis results, which can be achieved through a common identifier (such as a sub-area ID); creating a new data structure (such as a database table or data frame) containing the connection information in the topological relationship network (such as adjacency, distance, etc.) and the task continuity analysis results (such as terrain similarity), ensuring that the relevant data of each sub-area is integrated. The criteria for evaluating the spatial connectivity of sub-areas include:

[0118] Count the number of adjacent sub-regions directly connected to each sub-region. The more the number, the better the spatial connectivity of the sub-region. Calculate the straight-line distance from the center point of the sub-region to the center points of all its adjacent sub-regions and find the average value. The shorter the average distance, the better the spatial connectivity. Calculate the clustering coefficient of the sub-region, that is, the probability that the adjacent sub-regions of the sub-region are also adjacent to each other. The higher the clustering coefficient, the stronger the spatial connectivity of the sub-region. The criteria for evaluating task continuity include:

[0119] Use terrain data (such as elevation, slope, etc.) to calculate the terrain differences between adjacent sub-regions. Set a similarity threshold. When the terrain difference is lower than the threshold, the two sub-regions are considered to be topographically continuous. Compare the distribution of land object types in adjacent sub-regions through remote sensing image classification or land use data. Similarity indicators (such as Jaccard similarity) can be used to measure the consistency of the land object characteristics of the two sub-regions. Use GIS network analysis tools to calculate the shortest path length or required time from the center of one sub-region to the center of another sub-region. The lower the accessibility index value, the better the accessibility between the two sub-regions and the higher the task continuity.

[0120] For each evaluation criterion, determine its relative importance according to the actual situation and assign corresponding weights to it. The weights can be determined by the Delphi method; since the numerical ranges and units of different evaluation criteria may be different, the data needs to be standardized so that they can be compared and weighted summed on the same scale; for each sub-area, the weighted summation method is used to calculate the comprehensive evaluation value according to the defined evaluation criteria and corresponding weights. The specific calculation method is to multiply the score of each criterion by its weight, and then add the weighted scores of all criteria to obtain the comprehensive evaluation value of the sub-area; according to the size of the comprehensive evaluation value, the sub-areas are sorted and classified, and different thresholds or levels can be set according to actual needs to better understand and apply the evaluation results.

[0121] According to the comprehensive evaluation results, a grouping strategy is formulated to obtain relevant sub-region groups, including: determining the number of groups to be divided into according to actual needs and resource constraints; using a clustering algorithm (such as K-means) to divide the sub-regions into a specified number of groups according to the comprehensive evaluation results, ensuring that the sub-regions in the same group have high spatial and task correlation; optimizing and adjusting the preliminary grouping results, considering the constraints of the actual collection task (such as time, cost, staffing, etc.), to ensure the feasibility and effectiveness of the grouping strategy; outputting the final grouping results in an appropriate format, including a list of sub-regions in each group and related spatial visualization representations (such as group annotations on a map).

[0122] In an embodiment of the present invention, by identifying and grouping sub-areas that are adjacent or close in space, the collection path can be optimized, the moving time and cost of the agent between different sub-areas can be reduced, and the overall collection efficiency can be improved. Evaluating and considering the continuity of task execution between different sub-areas helps to ensure the smooth transition of collection tasks between adjacent sub-areas, reduce interruptions and duplication of work, and improve the consistency and efficiency of task execution. Grouping sub-areas based on comprehensive evaluation results can more reasonably allocate collection resources and agents, ensure that resources are effectively utilized between sub-areas with similar characteristics or needs, and avoid resource waste. The grouping strategy provides a flexible decision-making framework, which can dynamically adjust the grouping method of sub-areas according to different collection goals, task requirements and resource conditions to adapt to changing environments and needs. By comprehensively considering the spatial correlation and task continuity of sub-areas, data collection routes and strategies can be planned more specifically, so as to obtain more comprehensive, accurate and representative data, and improve data quality and the reliability of analysis results. For collection tasks that require cross-sub-area collaboration, the grouping strategy helps to promote collaborative cooperation between different sub-areas, realize resource sharing and complementary advantages, and promote balanced development of the entire region.

[0123] In a preferred embodiment of the present invention, polarization data of the target terrain area is extracted from the remote sensing data; and the polarization degree of the terrain area surface is calculated by analyzing the polarization data, including:

[0124] Obtaining raw data containing polarization information from remote sensing satellites or ground-based remote sensing equipment;

[0125] Filter out the data set covering the target terrain area from the original data, specifically including: using remote sensing data processing software to cut out the data subset covering the area from the original data according to the coordinate range of the target area; pre-processing the filtered data, such as denoising, radiation calibration, etc., to ensure data quality;

[0126] According to the characteristics of remote sensing data, the polarization band containing polarization information is determined, specifically including: establishing a connection with the remote sensing data storage system to obtain the metadata or related documents of the data; reading and parsing the metadata file to extract key information about the data band setting; in the parsed metadata, scanning the description information of all bands, including key parameters such as band name, central wavelength, bandwidth, etc.; recording the characteristics of each band and constructing a band information database or data table; according to the principle of polarization measurement, setting the identification standard of the polarization band, such as a specific band name mode, central wavelength range or bandwidth threshold, using the above standards, querying and matching in the band information database, and preliminarily screening out the bands that may contain polarization information.

[0127] For each polarization band, extract the corresponding polarization component, and combine the polarization bands into a polarization data set. For each area, calculate its polarization degree, including: receiving instructions for processing remote sensing data, including the location, format and processing requirements of the data; reading the original remote sensing data from the specified location according to the instructions, parsing the data structure, and confirming the data integrity and readability; scanning all bands in the original data, identifying the bands containing polarization information through built-in parameters and standards, and marking the identified polarization bands; configuring the initial parameters required for Stokes parameter extraction, such as polarization angle, according to the characteristics of remote sensing data and polarization bands degree resolution, noise threshold, etc.; perform radiometric calibration and atmospheric correction on the marked polarization band to eliminate the influence of sensor response differences and atmospheric scattering on polarization information; apply image processing techniques, such as filtering or smoothing, to reduce random noise in the data and improve the accuracy of polarization component extraction; construct the preprocessed polarization band data into a data matrix according to the requirements of the Stokes parameter extraction method, and each pixel corresponds to a data vector; according to the definition and mathematical formula of the Stokes parameter, operate on each data vector in the data matrix to calculate the four components of the Stokes parameter (I, Q, U, V). Among them, I represents the total light intensity, Q and U represent the linear polarization components in the horizontal and vertical directions respectively, and V represents the circular polarization component; from the calculated Stokes parameters, extract the horizontal polarization component (related to the Q parameter) and the vertical polarization component (related to the U parameter), and further convert the Q and U components into horizontal and vertical polarization intensities according to the polarization angle.

[0128] The extracted polarization components are integrated in the order of the original bands; the integrated data is converted into an internal processing format to optimize storage and computing efficiency; a new data structure is created to store the integrated and converted polarization components to form a polarization data set; for each region, the degree of polarization is calculated using the polarization degree formula based on the value of the polarization component, which can help understand the polarization characteristics of different regions.

[0129] Assign the degree of polarization to the corresponding area and generate a degree of polarization map, which specifically includes determining the analysis unit (such as pixel, geographic grid, etc.), which is the basic unit for assigning the degree of polarization value; associate the degree of polarization value of each area with its corresponding analysis unit, which can be achieved by creating a data table, using a spatial database or establishing an index, etc., to ensure that each analysis unit can accurately correspond to a degree of polarization value; start and initialize the GIS software or remote sensing image processing software; import the associated degree of polarization data into the image processing software; formulate a suitable color mapping scheme based on the value range and distribution characteristics of the degree of polarization. This involves selecting an appropriate color space (such as RGB, HSV, etc.) and setting the color or grayscale level corresponding to different degree of polarization values; create a new layer in the image processing software to display the degree of polarization map. According to the color mapping scheme, the degree of polarization value of each analysis unit is converted to the corresponding color or grayscale level and rendered on the layer; adjust the image brightness, contrast, resolution and other parameters as needed to optimize the visual effect and readability of the degree of polarization map; after the degree of polarization map is generated, preview and check the integrity and accuracy of the image. Convert polarization images to commonly used image formats (such as JPEG, PNG, TIFF, etc.) and save them to a specified storage location. At the same time, they can also be saved in software-specific data formats for further analysis and processing in the same software environment.

[0130] In an embodiment of the present invention, polarization data can provide additional information about the characteristics of surface reflected light, which helps to more accurately identify and distinguish different terrain features (such as water bodies, vegetation, bare soil, buildings, etc.). Through the calculation of polarization degree, these features can be further quantified, thereby improving the accuracy of terrain classification and identification. Traditional remote sensing images mainly rely on spectral information for interpretation, and the addition of polarization information provides a new dimension for image interpretation. The polarization degree map can be used as an auxiliary tool to help interpreters more intuitively understand the surface cover type and material properties, and reduce the uncertainty and subjectivity in the interpretation process. Polarization data is highly sensitive to environmental changes (such as soil moisture, vegetation health, etc.). By regularly acquiring and analyzing polarization data, dynamic monitoring of environmental parameters can be achieved.

[0131] In a preferred embodiment of the present invention, the incident angle and azimuth of the incident light are estimated according to the shooting angle and the sun position of the remote sensing data; a first mapping relationship between an initial normal vector and the incident angle and azimuth of the incident light is constructed using the information of the degree of polarization, the incident angle and the azimuth; the zenith angle and azimuth of the sun are calculated according to the shooting time and the geographical location information of the remote sensing data; and a second mapping relationship between a light source vector and the zenith angle and the azimuth is constructed using the zenith angle and the azimuth of the sun, including:

[0132] According to the format of remote sensing data (such as GeoTIFF, NetCDF, etc.), use the corresponding library (such as GDAL, NetCDF4, etc.) to read metadata; find information related to the sensor shooting angle in the metadata, such as "pitch", "roll", "yaw", etc., which describe the posture of the sensor at the time of shooting; store these angle information in variables. Extract the shooting time (including date and time) and geographical location information (longitude and latitude) of remote sensing data from metadata; use astronomical algorithms (such as SPA algorithm), input the shooting time and geographical location, and obtain the altitude and azimuth of the sun at that moment; according to the geographical location and shooting angle of remote sensing data, establish a local coordinate system, in which the Z axis is perpendicular to the ground, and the X axis and Y axis are on the ground plane; use the altitude and azimuth of the sun to calculate the three-dimensional vector of the sun's rays in the local coordinate system; calculate the three-dimensional vector of the sensor's line of sight according to the pitch, roll and yaw angles of the sensor; use the dot product formula of three-dimensional space vectors to calculate the angle between the sun's ray vector and the sensor's line of sight vector, which is the angle of incidence. Project vectors to the ground plane: Project the sun's ray vector and the sensor's sight line vector onto the XY plane (i.e., the ground plane); use the dot product and cross product formulas of two-dimensional vectors to calculate the angle between the two projected vectors, which is the azimuth.

[0133] According to the previous steps, the data sets of polarization degree, incident angle and azimuth angle are calculated; check whether there are missing values, outliers or erroneous data in the data set. For missing values, you can consider using interpolation or other methods to fill them; for outliers, you need to verify and correct them; among them, each polarization degree value is matched with the corresponding incident angle and azimuth angle values, which means that each data point should contain information about these three parameters; organize the data into a format suitable for subsequent processing, for example, you can create a data table or array containing three columns (polarization degree, incident angle, azimuth angle).

[0134] According to the known medium properties and incident angle, the refractive index is calculated using the Fresnel equation; the refractive index and incident angle are used to further calculate the reflection coefficient, which describes the reflection behavior of light on the medium interface; the calculated intermediate variables (such as refractive index, reflection coefficient) and input parameters (polarization degree, incident angle, azimuth) are substituted into the Fresnel equation; the components of the initial normal vector corresponding to the input parameters are obtained by solving the equation; the calculated initial normal vector components are combined into a three-dimensional vector. When it comes to specific applications, the following process can be used:

[0135] The refractive index is calculated based on the properties of the two media (such as optical density or the speed of light in the media). For a given angle of incidence and medium properties, the refractive index can be calculated using Snell's law (also called the law of refraction). The formula for Snell's law is:

[0136] ,in, and are the refractive indices of the two media, is the angle of incidence, is the angle of refraction, from which we can find the angle of refraction The calculation formula is:

[0137] ;

[0138] Calculate the reflection coefficient, where the reflection coefficient describes the proportion of energy reflected at the interface of two media; use the Fresnel equation to calculate the reflection coefficient. For unpolarized light, the Fresnel equation gives the relationship between the reflection coefficient and the angle of incidence and the refractive index of the two media; the reflection coefficient R can be calculated using the Fresnel equation. The (transverse electric) polarization calculation formula is:

[0139] ;

[0140] for The (transverse magnetic) polarization calculation formula is:

[0141] ;

[0142] Calculate the components of the initial normal vector. The normal vector is perpendicular to the incident medium interface. For a planar interface, the normal vector n It can be expressed as:

[0143] ;

[0144] In the two-dimensional case, if the angle of incidence is xy In the plane, then:

[0145] n x =0; n y =0; n z =1 (pointing out of the incident plane);

[0146] For the three-dimensional case, assume that the incident plane is xz The planes coincide and the normal vector can be written as:

[0147] ;

[0148] Combine the components calculated above to form the final normal vector N, where N = ( n x , n y , n z );in, ; ; .

[0149] The shooting time and geographic location information are extracted from the metadata of remote sensing data, and the zenith angle and azimuth angle of the sun at that moment are calculated based on the shooting time and geographic location using astronomical algorithms. The zenith angle is the angle between the sun's rays and the vertical line on the ground, and the azimuth angle is the angle between the projection of the sun's rays on the ground and the true north direction.

[0150] According to the sun's zenith angle θ and azimuth , we can define a three-dimensional vector to represent the direction of the sun's light source. In the Cartesian coordinate system, this three-dimensional vector S can be expressed as: S =( S x , S y , S z );in:

[0151] ;

[0152] ;

[0153] ;

[0154] Here, θ is the zenith angle, which represents the angle between the sun’s rays and the vertical line on the ground, ranging from 0° to 180°; The azimuth angle represents the angle between the projection of the sun's rays on the ground and the true north direction, ranging from 0° to 360°. Determined, so from S to (θ, ) is a direct mapping. This mapping can be achieved by the following formula:

[0155] θ=arccos( S z );

[0156] ;

[0157] Apply the constructed mapping relationship to remote sensing data. For each set of remote sensing data, extract its shooting time and geographical location information, and then use astronomical algorithms to calculate the corresponding zenith angle and azimuth angle. Then, use these angle values ​​to construct the light source vector S, and reversely calculate the new zenith angle and azimuth angle through the above mapping relationship; S x Indicated in x Components on the axis; S y Indicated in y Components on the axis; S z Indicated in z Components on the axis.

[0158] In a preferred embodiment of the present invention, the reflected light intensity of each area is calculated according to the normal vector information and the light source vector in the first mapping relationship, including:

[0159] The normal vector information of each area is obtained from the first mapping relationship, and the light source vector is calculated according to the zenith angle and azimuth angle of the sun, which specifically includes: traversing each area in the first mapping relationship and extracting the corresponding normal vector data; ensuring that the normal vector is a unit vector, that is, its length (or modulus) is 1. If not, it needs to be normalized; the zenith angle is the angle between the sun's rays and the ground normal, and the azimuth angle is the angle between the projection of the sun's rays on the horizontal plane and a certain reference direction (such as due north). These angle values ​​can be obtained through real-time weather data, time and date calculations, or user input; based on the sun's zenith angle and azimuth angle, a light source vector representing the direction of the sun's rays can be constructed, and the zenith angle and azimuth angle can be converted into coordinates in three-dimensional space. For example, in a right-handed coordinate system assuming that the Y axis is in the upward direction, the X axis is in the due east direction, and the Z axis is in the due north direction:

[0160] The X component of the light source vector = sin(zenith angle) × cos(azimuth angle); the Y component of the light source vector = -cos(zenith angle); the Z component of the light source vector = sin(zenith angle) × sin(azimuth angle);

[0161] According to the normal vector information in the first mapping relationship and the calculated light source vector, Calculate the reflected light intensity for each area ;

[0162] in, Indicates the direction of incident light; Indicates the viewing direction; represents the normal vector; represents the roughness parameter; is the half vector of incident and reflected light; where, ; is a parameter related to surface roughness and its value is ;in, ;in, is the reflectivity at normal incidence.

[0163] In the embodiment of the present invention, by combining the normal vector information and the light source vector, the calculation method can accurately simulate the reflected light intensity of different areas under different lighting conditions. The calculation method achieves efficient calculation performance by optimizing the mathematical formula and algorithm structure. This means that when processing a large number of areas or complex scenes, a relatively fast calculation speed can still be maintained.

[0164] In a preferred embodiment of the present invention, based on the calculation of the reflected light intensity and remote sensing data of each area, the prior normal vector information of the surface of the target terrain area is obtained through inversion calculation, including:

[0165] Obtain the surface reflectivity data of the target terrain area, and calculate the relationship between the incident angle and the reflection angle through the sun position, sensor viewing angle and surface geometric relationship, specifically including: obtaining the surface reflectivity data of the target terrain area from remote sensing satellites or ground measurement equipment, these data are represented in the form of images or numerical matrices, and each pixel or grid point corresponds to a surface reflectivity value; preprocessing the acquired reflectivity data, including denoising, correction and standardization operations; storing the processed reflectivity data as an appropriate data structure; obtaining or calculating the sun position information, including the zenith angle and azimuth angle, which describe the direction of the sun's rays relative to the surface; obtaining the sensor's viewing angle information, that is, the direction and angle when the sensor observes the surface; for each point or pixel in the target terrain area, according to its surface geometric features (such as elevation, slope, slope aspect, etc.), combined with the sun position and sensor viewing angle, the geometric relationship is used to calculate the incident angle and reflection angle of the point.

[0166] The residual is defined as the difference between the observed reflected light intensity and the light intensity predicted by the parameterized reflection process, which includes: reading the reflected light intensity value of the point or pixel from remote sensing data or ground measurement equipment. This value is usually a digital brightness or radiation value, representing the intensity of light reflected from the point; storing these observations in a data structure, such as an array or list; selecting a physically based illumination model, such as the Lambertian model, which can predict the reflected light intensity based on the surface material properties and lighting conditions; the reflected light intensity prediction formula of the Lambertian model is: H = k d ×W l ×cos(θ i ), where H is the reflected light intensity, k d is the diffuse reflectance of the surface, W l is the intensity of the incident light, θi It is the angle between the incident light and the surface normal vector. For each observation point or pixel, an initial normal vector can be assumed. For example, (0, 0, 1) indicates that the surface is horizontal at the beginning. In the subsequent optimization process, this normal vector will be adjusted to better fit the observed data.

[0167] Diffuse reflectance k d This coefficient describes the diffuse reflection ability of the surface, and an initial value can be set according to the material type of the surface. For example, for darker surfaces such as sand, a lower initial value can be set; for brighter surfaces such as snow, a higher initial value can be set.

[0168] Surface roughness (if the model takes this into account). Surface roughness affects how light is scattered. An initial roughness value can be set based on previous surface measurements or experience.

[0169] Based on the position of the sun (zenith angle and azimuth) and the surface geometry of the observation point or pixel (such as elevation, longitude and latitude, etc.), the incident direction of the sun's rays relative to the surface can be calculated. This usually involves converting the sun's position from an astronomical coordinate system to a local surface coordinate system.

[0170] Reflected light direction: In the Lambertian model, reflected light is evenly scattered. In order to match the sensor's viewing angle, the direction from the observation point or pixel to the sensor can be calculated and regarded as the observation direction of the reflected light. The Lambertian model and the calculated illumination direction are used in combination with the reflectivity of the surface and other material properties to predict the reflected light intensity at that point; the observed reflected light intensity value is compared with the reflected light intensity value predicted by the Lambertian model, and the difference between the two is calculated, that is, the residual, which represents the difference between the model prediction and the actual observation.

[0171] Construct the objective function of the residual sum of squares, set the initial values ​​for the unknown parameters in the reflection process, and the objective function The calculation formula is:

[0172] ;

[0173] in, is the number of observation points; At the observation point The observed intensity of reflected light; is the weight of each observation point; It is an observation point The standard deviation of the observation error; is the regularization parameter; is the angle between the incident light and the surface normal; is the azimuth angle between the incident light and the reflected light; Represents the index value; is the reflectivity of the surface; the gradient descent method is used to minimize the residual sum of squares through iterative calculation, and the parameter value is adjusted in each iteration, specifically including: using the gradient descent method as an optimization algorithm, and minimizing the objective function of the residual sum of squares through iterative calculation; in each iteration, the gradient of the objective function with respect to the unknown parameter is calculated (that is, the direction in which the function value changes fastest); according to the calculated gradient information, the value of the unknown parameter is adjusted so that the objective function moves toward the minimum value, and the iterative calculation is repeated until the stop condition is met (such as reaching a preset number of iterations, the objective function value converges, etc.); the residual sum of squares between the predicted data and the observed data after each iteration is calculated, and when the change in the residual sum of squares reaches a preset number of iterations, the iteration is stopped, and the parameter value in the current iteration step is output as the final solution, wherein the final solution of the normal vector is the prior normal vector information of the surface of the target terrain area; when applied specifically, the following steps are specifically included:

[0174] Set the number of observation points N; for each observation point i, initialize the observed reflected light intensity , weight wi, observation error standard deviation σi; initialize regularization parameter λ; set surface reflectivity An initial value for the unknown surface normal vector, which is usually a vector perpendicular to the surface, such as (0, 0, 1).

[0175] Calculate the direction of incident light based on the sun position and the surface geometry of the observation point; calculate the direction of reflected light based on the sensor position and the surface geometry of the observation point; calculate the angle between the incident light and the surface normal from the directions of the incident light and reflected light and the azimuth angle between the incident light and the reflected light .

[0176] The Lambertian reflectance model is used to predict the reflected light intensity Ii at each observation point, which involves using the surface reflectance , Angle , azimuth .

[0177] For each observation point i, calculate the observed reflected light intensity Reflected light intensity predicted by the model The difference between the two is the residual; the residual is divided by the standard deviation of the observed error , get the standardized residual; square the standardized residual and multiply it by the corresponding weight .

[0178] The weighted standardized residual squares of all observation points are summed to obtain the value of the first term; the value of the regularization term is calculated, that is, λ multiplied by the expression consisting of the cosine value of the surface reflectivity, angle, and azimuth; the first term and the regularization term are added to obtain the objective function The value of; use the gradient descent method to minimize the objective function ; In each iteration, calculate the gradient of the objective function with respect to unknown parameters (such as surface normal vector, reflectivity, etc.); update the value of the unknown parameter according to the gradient information to reduce the value of the objective function; repeat the iterative process until the preset stopping condition is reached (such as the number of iterations, the change in the objective function value is less than the threshold, etc.); when the optimization process stops, output the unknown parameter values ​​in the current iteration step as the final solution. These parameter values ​​include the optimized surface normal vector, reflectivity, etc., which can best fit the observed data and meet the physical constraints.

[0179] In an embodiment of the present invention, the surface normal vector information of the target terrain area can be obtained more accurately through inversion calculation. The acquisition of prior normal vector information provides more valuable auxiliary information for the interpretation of remote sensing data. Combining the reflected light intensity and remote sensing data, we can have a deeper understanding of the nature, state and distribution of surface materials, thereby improving the accuracy and efficiency of remote sensing data interpretation. In the illumination model, the normal vector is one of the key parameters. The first term in the objective function weights the square of the residual of each observation point by weight and the standard deviation of the observation error, so as to more accurately reflect the uncertainty of the observation data. At the same time, the second term introduces physical constraints (such as the Lambertian reflection model) to ensure the physical rationality of the solution. By introducing the standard deviation of the observation error , the objective function can be more robust to noise or outliers during the optimization process. Points with large observation errors will receive smaller weights in the optimization, thereby reducing their impact on the final solution. By setting reasonable initial values ​​for unknown parameters and combining optimization algorithms such as gradient descent, the minimum value of the objective function can be searched more effectively, improving the efficiency and stability of the optimization process.

[0180] In a preferred embodiment of the present invention, the normal vector information in the first mapping relationship is corrected based on the prior normal vector information to obtain the corrected normal vector information, including:

[0181] The normal vector information in the first mapping relationship of the target terrain area is obtained from the data storage, which is recorded as mapped_normals; at the same time, the prior normal vector information of the same area is obtained, which is recorded as prior_normals.

[0182] For each normal vector in mapped_normals, find the normal vector in prior_normals that corresponds to the spatial position; if the resolution of prior_normals is inconsistent with mapped_normals, make their resolution consistent through interpolation or resampling methods to ensure that each mapped normal vector has a corresponding prior normal vector.

[0183] Prior_normals is evaluated for quality, including checking the integrity, consistency, noise level, etc. of the data to obtain the evaluation results; according to the evaluation results, a quality score is assigned to each prior normal vector, denoted as quality_scores. The higher the quality score, the higher the credibility of the prior normal vector.

[0184] Combined with the source credibility of the prior normal vector, the correction weight of each prior normal vector is calculated, recorded as correction_weights; for each mapped normal vector and its corresponding prior normal vector, the difference vector difference_vector between them is calculated, that is:

[0185] difference_vector=prior_normal-mapped_normal.

[0186] Use the correction weight to scale the difference vector to get the correction value correction_vector, where:

[0187] correction_vector=correction_weight×difference_vector.

[0188] Add the correction correction_vector to the corresponding mapped normal vector to get the corrected normal vector corrected_normal, where:

[0189] corrected_normal=mapped_normal+correction_vector.

[0190] Combine all corrected normal vectors corrected_normal to form the corrected normal vector information.

[0191] In a preferred embodiment of the present invention, three-dimensional scene reconstruction of the target area is performed based on the corrected normal vector information and polarization data of different target terrain areas at a preset angle, including:

[0192] Load the corrected normal vector information and polarization data of different target terrain areas at preset angles from the data storage; create an empty 3D scene and set the appropriate coordinate system and spatial range; based on the corrected normal vector information and the preset shooting angle, use geometric relationships to convert the normal vector into point coordinates in 3D space, which constitute the preliminary 3D point cloud of the target terrain area; use the surface reflection information contained in the polarization data to make fine adjustments to the preliminary 3D point cloud; combine the constraints of the normal vector and polarization data to optimize the position and orientation of each point to more accurately reconstruct the geometry and details of the terrain surface.

[0193] Extract the grid structure of the terrain surface from the refined 3D point cloud. Use triangulation algorithms (such as Delaunay triangulation) to generate a continuous and seamless surface grid, map remote sensing images or other texture data onto the reconstructed terrain surface, and assign corresponding material maps to different areas according to the material properties of the terrain. Save the reconstructed 3D scene as a common 3D model file format (such as OBJ, FBX, etc.); use 3D visualization software or virtual reality equipment to display and interactively explore the reconstructed target terrain area.

[0194] The above is a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A three-dimensional terrain mapping system based on remote sensing technology, characterized in that: The system comprises: Data acquisition module, used to collect remote sensing data of different target terrain areas; The data processing module is used to extract polarization data of the target terrain area from the remote sensing data; calculate the polarization degree of the terrain area surface by analyzing the polarization data; estimate the incident angle and azimuth of the incident light according to the shooting angle of the remote sensing data and the position of the sun; use the information of the polarization degree, the incident angle and the azimuth to construct a first mapping relationship between an initial normal vector and the incident angle and azimuth of the incident light; calculate the zenith angle and azimuth of the sun according to the shooting time and geographical location information of the remote sensing data; use the zenith angle and azimuth of the sun to construct a second mapping relationship between a light source vector and the zenith angle and azimuth; calculate the reflected light intensity of each area according to the normal vector information and the light source vector in the first mapping relationship; obtain the prior normal vector information of the target terrain area surface through inversion calculation based on the calculation of the reflected light intensity of each area and the remote sensing data; correct the normal vector information in the first mapping relationship based on the prior normal vector information to obtain the corrected normal vector information; A three-dimensional reconstruction module is used to reconstruct the three-dimensional scene of the target area based on the corrected normal vector information and the polarization data of different target terrain areas at a preset angle; Extract polarization data of the target terrain area from remote sensing data; calculate the polarization degree of the terrain area surface by analyzing the polarization data, including: Obtaining raw data containing polarization information from remote sensing satellites or ground-based remote sensing equipment; Filtering the dataset covering the target terrain area from the original data; According to the characteristics of remote sensing data, the polarization band containing polarization information is determined; For each polarization band, the corresponding polarization component is extracted, and the polarization bands are combined into a polarization data set; For each region, calculate its degree of polarization; Assigning polarization degrees to corresponding regions to generate a polarization degree map; calculating the reflected light intensity of each region according to the normal vector information and the light source vector in the first mapping relationship, including: Obtain normal vector information of each region from the first mapping relationship, and calculate the light source vector according to the zenith angle and azimuth angle of the sun; According to the normal vector information in the first mapping relationship and the calculated light source vector, Calculate the reflected light intensity for each area ; in, Indicates the direction of incident light; Indicates the viewing direction; represents the normal vector; represents the roughness parameter; is the half vector of incident and reflected light; where, ; is a parameter related to surface roughness and its value is ;in, ;in, is the reflectivity at normal incidence; based on the calculation of the reflected light intensity and remote sensing data of each area, the prior normal vector information of the surface of the target terrain area is obtained through inversion calculation, including: Obtain the surface reflectance data of the target terrain area, and calculate the relationship between the incident angle and the reflection angle based on the sun position, sensor viewing angle and surface geometry; The residual is defined as the difference between the observed reflected light intensity and the light intensity predicted by the parameterized reflectance process; Construct the objective function of the residual sum of squares and set the initial values ​​for the unknown parameters in the reflection process; The gradient descent method is used to minimize the residual sum of squares through iterative calculation. In each iteration, the parameter value is adjusted; Calculate the residual sum of squares between the predicted data and the observed data after each iteration. When the change of the residual sum of squares reaches the preset number of iterations, stop the iteration and output the parameter value in the current iteration step as the final solution, where the final solution of the normal vector is the prior normal vector information of the surface of the target terrain area; the objective function The calculation formula is: ； in, is the number of observation points; At the observation point The observed intensity of reflected light; is the weight of each observation point; It is an observation point The standard deviation of the observation error; is the regularization parameter; is the angle between the incident light and the surface normal; is the azimuth angle between the incident light and the reflected light; Represents the index value; is the reflectivity of the surface.

2. The three-dimensional terrain mapping system based on remote sensing technology according to claim 1 is characterized in that: Collect remote sensing data of different target terrain areas, including: Divide the target terrain area into multiple sub-areas, each of which is responsible for one or more agents; Assign initial remote sensing data collection parameters to each agent, including collection time, sensor configuration, and flight path; Each agent simulates remote sensing data collection in the responsible sub-area according to its current collection parameters and evaluates the performance of the collection parameters to obtain local evaluation results; Based on the information shared by all agents, a global estimation is performed to evaluate the overall performance of the current acquisition parameters to obtain a global estimation result; According to the global estimation results, the acquisition parameters of each agent are updated, and the local estimation and global estimation are repeated until the preset number of iterations is reached to obtain an optimized remote sensing data acquisition plan; According to the optimized remote sensing data collection plan, remote sensing data of different target terrain areas are collected.

3. The three-dimensional terrain mapping system based on remote sensing technology according to claim 2 is characterized in that: The target terrain area is divided into multiple sub-areas, each of which is managed by one or more agents, including: Determine the scope and boundaries of the target terrain area; preliminarily divide the target terrain area into multiple sub-areas based on the terrain's height, slope, vegetation cover characteristics, and specific requirements for remote sensing data collection; For each initially divided sub-area, extract key terrain and ground features; Analyze key terrain and landform features, and classify sub-areas according to the similarity of key terrain and landform features to obtain classification results; According to the classification results, the spatial correlation between the various sub-regions and the continuity of the acquisition tasks are analyzed to obtain the sub-region groups with correlation; Each sub-region that has been classified and determined to be relevant is regarded as a node in the graph model, and the feasible paths between sub-regions are regarded as edges in the graph; Define the rules for ants to move in the graph, initialize the ant population, and each ant chooses the next sub-area to move to based on the pheromone concentration and heuristic information of the nodes around the current location. Ants release pheromones when selecting a path and update the pheromone concentration based on the length and quality of the path. When all ants have completed a path search, they will globally update the pheromone according to the length and quality of the overall path until the preset number of iterations is reached. After the iteration is completed, the optimized sub-area division scheme and the corresponding path planning results will be output; According to the optimized sub-area division scheme and the corresponding path planning results, the intelligent agent is assigned to the corresponding sub-area for remote sensing data collection.

4. The three-dimensional terrain mapping system based on remote sensing technology according to claim 3 is characterized in that: Analyze key terrain and landform features, and classify sub-areas according to the similarity of key terrain and landform features to obtain classification results, including: The terrain and object features of each sub-region are combined into a feature vector, which represents the comprehensive features of the sub-region; Determine the corresponding K value through cross-validation, prepare an empty result list for the sub-regions to be classified; for each sub-region to be classified, perform the following steps: Calculate the distance between the feature vector of the current sub-region and the feature vectors of all sub-regions in the training set; Sort the distances and select the corresponding K sub-regions as the neighbors of the current sub-region to obtain the classification result.

5. The three-dimensional terrain mapping system based on remote sensing technology according to claim 4 is characterized in that: According to the classification results, the spatial correlation between the various sub-regions and the continuity of the acquisition tasks are analyzed to obtain relevant sub-region groups, including: Analyze the adjacency relationship between each sub-region, and for sub-regions that are not directly adjacent, calculate the spatial distance between the sub-regions; Based on the adjacency relationship and spatial distance, a topological relationship network between sub-regions is established, and the topological relationship network represents the spatial connectivity of sub-regions; According to the goal of the collection task, evaluate the continuity of task execution between different sub-areas to obtain the analysis results of the continuity of the collection task; Based on the analysis results of the topological relationship network and the continuity of the acquisition task, a comprehensive evaluation is conducted on the sub-areas to obtain a comprehensive evaluation result; Based on the comprehensive evaluation results, a grouping strategy is formulated to obtain relevant sub-region groups.

6. A three-dimensional terrain mapping method based on remote sensing technology, characterized in that: Applicable to the system according to any one of claims 1 to 5, comprising: Collect remote sensing data of different target terrain areas; Extract polarization data of the target terrain area from remote sensing data; calculate the polarization degree of the terrain area surface by analyzing the polarization data; estimate the incident angle and azimuth of the incident light according to the shooting angle of the remote sensing data and the position of the sun; use the information of the polarization degree, incident angle and azimuth to construct a first mapping relationship between an initial normal vector and the incident angle and azimuth of the incident light; calculate the zenith angle and azimuth of the sun according to the shooting time and geographical location information of the remote sensing data; use the zenith angle and azimuth of the sun to construct a second mapping relationship between a light source vector and the zenith angle and azimuth; calculate the reflected light intensity of each area according to the normal vector information and the light source vector in the first mapping relationship; obtain the prior normal vector information of the surface of the target terrain area through inversion calculation based on the calculation of the reflected light intensity of each area and the remote sensing data; correct the normal vector information in the first mapping relationship based on the prior normal vector information to obtain the corrected normal vector information; The three-dimensional scene of the target area is reconstructed based on the corrected normal vector information and the polarization data of different target terrain areas at preset angles.

Images