Automatic surveying and mapping analysis method and system based on remote sensing data of unmanned aerial vehicle

Abstract

The invention belongs to the technical field of geospatial information surveying and mapping, and particularly relates to an automatic surveying and mapping analysis method and system based on unmanned aerial vehicle remote sensing data. The method comprises the following steps: firstly, acquiring multi-stage unmanned aerial vehicle remote sensing data, optimizing a strong backlight scene by using a multi-exposure fusion algorithm, correcting reflectivity in a rainstorm ponding area by using a dark target method, denoising laser point cloud by using neighborhood statistical filtering, and screening a high-definition surveying and mapping area with a complete edge; secondly, associating time sequence information of the same spatial position, calculating a vegetation quarterly growth rate through NDVI index time sequence fitting, calculating an earth surface accumulated settlement amount and an annual settlement rate through a radar interference measurement technology, and constructing a dynamic change feature set; and finally, a threshold value is preset to verify the characteristic parameters, unmatched artificial buildings are identified, and an automatic report is generated after manual intervention calibration is triggered. According to the method, scene-divided data optimization is innovatively realized to break through the complex scene surveying and mapping precision bottleneck, and meanwhile, the problem of unbalance between automation and flexibility is solved by constructing an intelligent verification and manual intervention cooperation mechanism.

Description

Technical Field

[0001] This invention relates to the field of geospatial information mapping technology, specifically to an automated mapping and analysis method and system based on UAV remote sensing data. Background Technology

[0002] With the rapid advancement of the social economy, the demand for geospatial information in various industries has shifted from traditional static data acquisition to dynamic monitoring, real-time analysis, and precise feedback. This necessitates the ability to quickly capture the changing characteristics of spatial elements and accurately extract the underlying information from the data, providing timely and effective data support for various scenarios. The rise and widespread adoption of UAV remote sensing technology, with its unique advantages of flexibility, wide coverage, and high data acquisition efficiency, has become an important means of geospatial information collection, providing strong support for meeting the information needs of multiple scenarios and high frequencies. Constructing an efficient and accurate automated analysis system is key to fully releasing the value of UAV data and adapting it to the modern information application needs of various industries.

[0003] However, while existing UAV remote sensing automated mapping systems perform reasonably well in conventional scenarios, their analysis accuracy drops significantly when faced with challenging scenarios such as strong backlighting, flooded areas after heavy rain, and complex terrain at high altitudes. For example, in strong backlighting environments, if optical remote sensing images are overexposed or underexposed, the system cannot effectively extract features or misidentify ground features. Due to the lack of a targeted adaptive adjustment mechanism, the mapping results are unreliable in complex scenarios. Secondly, existing UAV remote sensing mapping automated analysis either relies excessively on preset programs or requires frequent manual intervention to adjust parameters. For instance, when an unknown large man-made structure suddenly appears in the mapping area, the automated system cannot accurately identify and classify it, or manual modification of the analysis rules is required, affecting the continuity of mapping.

[0004] Therefore, this invention proposes an automated mapping and analysis method and system based on UAV remote sensing data. Summary of the Invention

[0005] To address the technical problems mentioned in the background art, such as the lack of adaptive robustness to edge scenes and poor integration between automated decision-making and manual intervention in existing UAV remote sensing automated mapping systems, the purpose of this invention is to provide an automated mapping analysis method and system based on UAV remote sensing data.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: Automated mapping and analysis methods based on UAV remote sensing data include: S1: Acquire multiple periods of UAV remote sensing data and perform scenario-based UAV remote sensing data optimization processing, extract the mapping areas that meet the conditions in each period of data, and output the multi-period remote sensing dataset; S2: Based on the aforementioned multi-period remote sensing dataset, associate the temporal characteristics of the same spatial location in the surveyed area, and through... NDVI Exponential time series fitting is used to calculate vegetation growth rate, radar interferometry is used to calculate land subsidence and subsidence trend, and a dynamic change feature set is output. S3: Based on the dynamic change feature set, preset the normal fluctuation range of vegetation growth rate and the safety threshold of land subsidence and verify them to determine whether it is necessary to trigger the manual interface. After manual confirmation, update the feature threshold or supplement the building feature library and output an automated surveying and analysis report.

[0007] Furthermore, the data obtained from multiple phases of UAV remote sensing are as follows: S101: The acquired multi-period UAV remote sensing data set is... ;in, This represents the total number of data collection periods, and ; For the first The first phase of the UAV remote sensing dataset includes the first phase. A subset of optical images from the period LiDAR point cloud subset and UAV remote sensing data acquisition timestamp .

[0008] Furthermore, the optimization and processing flow for drone remote sensing data in different scenarios is as follows: S102: For optical images with strong backlighting, an optimized optical image is obtained using a multi-exposure fusion algorithm, with the following formula: in, These are the planar coordinates of the pixel. This represents the total number of optical image frames. For optical image frame indexing; To optimize the optical image at the pixel level The grayscale value at that location; For the first Issue No. Frame optical images in The original grayscale value of a pixel; For the first Frame optical images in The pixel weights satisfy the conditions ; S103: To correct the reflectivity deviation in flooded areas after heavy rain, the dark target method is used. The formula is: in, To correct the water body at wavelength Reflectance at that location; To correct the water body at wavelength The original reflectivity at that location is determined by of DN Value conversion, and , wavelength Image at the location DN value, wavelength Saturated radiance at that location For corresponding saturation DN value; For dark targets at wavelength Reflectance at that location; S104: For isolated noise points in the lidar point cloud, the processing steps are as follows: S1041: Calculate the subset of the lidar point cloud The Middle Points The set of neighborhood points is given by the formula: in, for Surrounding distance ≤ The set of all neighboring points; The neighborhood search radius; for The Middle One point; For point 3D coordinates; S1042: Based on the neighborhood point set Calculation points Average distance to all points in the neighborhood To reflect the point The degree of fit with the surrounding terrain was determined, and the average distance of all points in the neighborhood was calculated. mean and standard deviation ; S1043: Based on the average distance of all points in the neighborhood mean and standard deviation Noise point identification and removal: If Judgment point These are noise points, and they should be removed; if Judgment point Points that are considered valid are retained; the final retained valid points constitute a point cloud. The terrain fitting vertical error must be satisfied. The formula is: in, This represents the number of points in the filtered point cloud. For point Vertical coordinates; This represents the actual terrain elevation. S1044: The optimized optical image data is used to screen the mapping area, which must simultaneously meet two quantitative indicators: sharpness and feature edge integrity. The clarity is measured by a structural similarity index. SSIM Assessment, requirements , Provide a clear reference image for the same area; the integrity of the ground feature edges is determined by edge matching degree. M Assessment, requirements ; Finally, the data of the surveyed area are integrated into... ,in, , corresponding to the number The optimized optical image, filtered lidar point cloud, UAV remote sensing data acquisition timestamp, optical image structure similarity index, UAV remote sensing data edge matching degree, and terrain fitting vertical error are analyzed.

[0009] Furthermore, for any spatial coordinate point within the surveying area... Perform the following matching operation: S201: From the optimized optical image Extract coordinates Corresponding pixel grayscale value The focus is on filtering the grayscale values ​​in the near-infrared and red light bands, and then obtaining the coordinates from the filtered lidar point cloud. Mean vertical coordinates of the surrounding effective point cloud , as coordinates No. The average topographic elevation over a period of time; to obtain coordinates. In the Time series datasets in the data : The process for calculating the vegetation growth rate is as follows: S202: Calculate the grayscale values ​​of the near-infrared and red bands of the optimized optical image. Period coordinates of NDVI The value, the formula is: in, For the first Period coordinates Optimized grayscale values ​​in the near-infrared band; For the first Period coordinates Optimized grayscale values ​​in the red band; For the first Period coordinates The normalized vegetation index at the location, with a value range of 100%. Negative values ​​indicate non-vegetated areas, while larger positive values ​​indicate more vigorous vegetation growth. Then, time difference conversion was performed, and the collection timestamps for each period were formatted as year, month, and day. Converted to time difference based on period 1 The unit is quarters, based on the time difference mentioned above. The independent variable is the normalized vegetation index. Using the dependent variable, construct a linear fitting model: in, For any time difference Below NDVI Predicted value; The slope of the fitted line, i.e., the coordinates Quarterly growth rate of vegetation at the location; The intercept of the fitted line represents the first period. of NDVI Initial value; Minimize the fitting error using the least squares method. Solve and And by calculating the goodness of fit To perform validity verification, the following requirements must be met. to indicate time sequence NDVI The data fits the linear model well, and the growth rate calculation results are reliable.

[0010] Furthermore, the process for calculating the surface subsidence and subsidence trend is as follows: S203: Based on the topographic elevation of Phase 1 Based on this, calculate the first... Period coordinates Cumulative settlement Negative values ​​represent settlement, and positive values ​​represent uplift. The formula is: in, For the first Period coordinates The average topographic elevation; This is the baseline elevation value for Phase 1; Time difference in quarters Further converted into time difference in years ; Subsequently For independent variable, Using the variable as the dependent variable, construct a linear fitting model and solve for the annual settlement rate. Negative values ​​represent average annual subsidence, and positive values ​​represent average annual uplift. The formula is: in, For any annual time difference, the cumulative settlement is predicted. The intercept of the fitted line represents the first period. The initial settlement; Similarly, the solution is obtained using the least squares method. and And by calculating the goodness of fit of the settlement To perform validity verification, the following requirements must be met. This indicates that the time-series settlement data fits the linear model well, and the calculation results of settlement rate and trend are reliable; The integration result is the dynamically changing feature set, and the process is as follows: S204: For all valid coordinate points within the surveyed area It must meet the following conditions simultaneously. , as well as Construct a dynamically changing feature set ,and , Represents the effective mapping area; where the feature set of a single coordinate point is... Defined as , This is the timestamp for the latest data collection.

[0011] Furthermore, regarding dynamically changing feature sets Each valid coordinate point The feature parameters are used to perform threshold checks on vegetation growth rate and land subsidence rate, respectively, to output the results of exceeding the limits. The process is as follows: S301: The vegetation growth rate threshold verification process is as follows: The vegetation growth rate threshold verification is performed using the actual quarterly vegetation growth rate at coordinate points. Relative deviation from the normal range Implementation; if Marked as , indicates that the limit has been exceeded; if Marked as , indicating normal; among them and All are labeled variables; S302: The surface subsidence rate threshold verification process is as follows: The surface subsidence rate threshold verification is performed by checking the coordinate points. The absolute value of the actual annual settlement rate is compared with the safety threshold for exceeding the limit. Calculation implementation; if Marked as , indicates that the limit has been exceeded; if Marked as This indicates normal; similarly, it includes... To mark variables; S303: For each coordinate point Define the comprehensive overlimit flag as Flag The formula is: in, Represents coordinate points At least one parameter is out of limit. This indicates that both parameters are normal; After completing the threshold verification of vegetation growth rate and land subsidence rate, the identification of unmatched artificial structures is carried out, and the process is as follows: S304: Based on the optimized optical image Extract coordinate points Surrounding building feature parameters Including the complexity of building outline Building height and spectral reflectance characteristics ; Subsequently, a comparison with the building feature database is performed, and a pre-set artificial building feature database is established. It includes the feature threshold range of common man-made structures, and ;in, The threshold range for building outline complexity. To minimize the contour complexity, The maximum contour complexity; For building height threshold range, Minimum building height ， Maximum building height; Minimum reflectivity Maximum reflectivity; Define the matching degree as ,like That is, at least two features do not match, and are marked as , indicates no matching artificial structure; otherwise, it is marked as , indicates that a match has been found; similarly, in Used to label variables.

[0012] Furthermore, when the coordinate point satisfy or When the manual intervention interface is triggered, the relevant parameters are updated after calibration, and the final verification result is output. The process is as follows: S305: The intervention trigger condition is all valid coordinate points within the survey area. or The point triggers a manual intervention interface with feature annotation, and the annotation content includes the point exceeding the limit. , Calculation results and unmatched buildings Details of the comparison between feature parameters and the feature library; then parameter calibration and updates are performed, including vegetation growth rate threshold calibration, land subsidence threshold calibration, and supplementation of the building feature library; S306: Based on the calibrated parameters, recalculate the comprehensive verification result for each coordinate point. The formula is: in, This indicates that it is normal; This indicates that the vegetation growth rate has exceeded the limit; This indicates that the surface subsidence has exceeded the limit; This indicates that both parameters have exceeded the limit; This indicates that no matching artificial structure was found. 、 , and All are recalculated marker values ​​after calibration; S307: Integrate all calculation results and verification information to generate an automated surveying and mapping analysis report. The automated surveying and mapping analysis report includes an overview of the surveying area, dynamic change statistics, verification result details, and recommended measures.

[0013] An automated mapping and analysis system based on UAV remote sensing data, wherein the system is applied to the automated mapping and analysis method based on UAV remote sensing data, the system includes the following steps: M1: Data acquisition module, which acquires remote sensing data from multiple periods through the UAV data acquisition unit and transmits it to the scene-specific optimization processing unit; it performs scene-specific processing such as multi-exposure fusion optimization, dark target reflectivity correction, and neighborhood analysis noise removal for strong backlight optical images, post-rainstorm water accumulation area images, and lidar point clouds respectively; then, through the region filtering unit, it extracts the mapping areas in each period of data that meet the requirements of sharpness ≥ 0.9 and ground feature edge integrity ≥ 0.85, and outputs multi-period remote sensing datasets; M2: Dynamic feature calculation module, which uses the time-series matching unit to associate time-series information of the same spatial location in the surveyed area from multiple remote sensing datasets to construct a coordinate point time-series dataset; the vegetation rate calculation unit is based on... NDVIThe exponential time series fitting model solves and outputs the quarterly vegetation growth rate; the subsidence trend calculation unit uses radar interferometry technology to calculate the cumulative surface subsidence and annual subsidence rate, integrates growth rate and subsidence trend data to form a dynamic change feature set and outputs it. M3: Feature verification and report generation module. The threshold verification unit presets the normal fluctuation range of vegetation growth rate and the safety threshold of land subsidence, and judges the parameters in the dynamic feature set that exceed the limit. The building identification unit extracts the building features around the coordinate point and compares them with the preset building feature library to identify unmatched artificial buildings. If there are buildings with parameters exceeding the limit or unmatched buildings, the manual intervention interface is triggered. After manual confirmation, the feature threshold is updated or the building feature library is supplemented. Finally, the report integration unit outputs an automated surveying and analysis report.

[0014] Compared with the prior art, the advantages of the present invention are as follows: 1. This invention adopts a scene-specific data optimization system, targeting three types of edge scenes: strong backlight, rain-affected areas, and high-altitude complex terrain. It deploys multi-exposure fusion algorithms, dark target reflectivity correction models, and laser point cloud neighborhood statistical filtering technology, respectively. This solves the core problems of overexposure and underexposure of optical images, spectral distortion in water-affected areas, and noise in point clouds of complex terrain. This ensures that the system maintains the same high reliability in various extreme scenarios as in conventional scenarios, and comprehensively improves the adaptability of the surveying and mapping system to complex environments.

[0015] 2. This invention designs an intelligent verification and dynamic rule update module. When special scenarios such as unknown man-made buildings are identified, the system automatically triggers the feature verification process. Through the manual intervention interface, parameter calibration and dynamic updates of analysis rules are realized. This not only avoids the drawbacks of preset program jams, but also greatly reduces the frequency and cost of manual intervention. It achieves the organic unity of the continuity of automated analysis and the accuracy of manual intervention, and achieves the goal of fully automatic and low-intervention efficient surveying.

[0016] 3. This invention utilizes time-series analysis techniques to integrate multi-period remote sensing data, and constructs coordinate points... NDVI The time-series fitting model and the surface subsidence interferometry algorithm can accurately extract dynamic change features such as vegetation growth rate and cumulative surface subsidence. Combined with an automated report generation system, it realizes intelligent processing of the entire process from data collection and feature extraction to analysis report output. It provides dynamic decision-making basis with both timeliness and accuracy for fields such as terrain monitoring, ecological management, and urban planning, and promotes the technological leap of surveying and mapping analysis from static data acquisition to dynamic trend prediction. Attached Figure Description

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

[0018] Picture 1 This is a schematic diagram of the workflow of the method of the present invention; Picture 2 This is a schematic diagram of the multi-period remote sensing dataset acquisition process of the present invention; Picture 3 This is a schematic diagram of the automated surveying and analysis report acquisition process of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0020] To achieve the above objectives, the present invention provides an automated mapping and analysis method based on UAV remote sensing data, such as... Picture 1-Picture 3 As shown, the method includes: S1: Acquire multiple periods of UAV remote sensing data and perform scenario-based UAV remote sensing data optimization processing, extract the mapping areas that meet the conditions in each period of data, and output the multi-period remote sensing dataset.

[0021] The following data was acquired from multiple UAV remote sensing datasets: S101: The acquired multi-period UAV remote sensing data set is... ;in, This represents the total number of data collection periods, and To ensure the effectiveness of time series analysis; For the first The first phase of the UAV remote sensing dataset includes the first phase. A subset of optical images from the period LiDAR point cloud subset and UAV remote sensing data acquisition timestamp ; No. The optical image subset mentioned in the period Represented as , For the first Total number of optical image frames in a period. For the first The first optical image subset Frame optical image; the first The lidar point cloud subset described in this period Represented as , For the first Total number of points in the lidar point cloud during the period For the first The first phase of the lidar point cloud sub-collection Each lidar point cloud, and , indicating the first Issue No. The three-dimensional coordinates of a LiDAR point cloud. Corresponding in sequence The east, north, and vertical coordinates of the point; All data is uniformly converted CGCS 2000 coordinate system, and synchronous data acquisition POS The data (position and attitude data) is used for lidar point cloud coordinate calibration; The workflow for optimizing and processing drone remote sensing data for different scenarios is as follows: S102: For the aforementioned subset of optical images The strong backlight optical image in the image is fused with multiple exposure frames to obtain the optimized optical image, as shown in the formula: in, These are the planar coordinates of the pixel. For optical image frame indexing; To optimize the optical image at the pixel level The grayscale value at that location; For the first Issue No. Frame optical images in The original grayscale value of a pixel; For the first Frame optical images in The weight of each pixel is calculated by combining image contrast, saturation, and brightness, and is determined based on the following conditions. ; Brightness standard deviation of the fused 8-bit format image To eliminate interference from excessively bright or dark backlighting; S103: To address the optical image reflectivity deviation caused by specular reflection in flooded areas after heavy rain, a dark target method is used for correction, with the following formula: in, To correct the water body at wavelength Reflectance at that location; To correct the water body at wavelength The original reflectivity at that location is determined by of DN Value conversion, and , wavelength Image at the location DN value, wavelength Saturated radiance at that location For corresponding saturation DN value; For dark targets (pure water in the shadow area) at wavelength Reflectance at (near-infrared band) ,at this time ); Corrected water accumulation area To avoid confusion between the reflectivity of water bodies and other land features; S104: Regarding the problem of isolated noise points generated in the lidar point cloud, the processing steps are as follows: S1041: Calculate the... The lidar point cloud subset described in this period The Middle Points The set of neighborhood points is given by the formula: in, for Surrounding distance ≤ The set of all neighboring points; The neighborhood search radius; for The Middle One point; For point The three-dimensional coordinates, corresponding in sequence The east, north, and vertical coordinates of the point; S1042: Based on the neighborhood point set Calculation points Average distance to all points in the neighborhood To reflect the point The formula for the degree of integration with the surrounding terrain is: in, For the set of neighborhood points The number of points contained therein, and To ensure the validity of statistics; Then the average distance of all points in the neighborhood is calculated. mean and standard deviation ; S1043: Based on the average distance of all points in the neighborhood mean and standard deviation Noise point identification and removal: If Judgment point These are noise points, and they should be removed; if Judgment point Points that are considered valid are retained; the final retained valid points constitute a point cloud. The terrain fitting vertical error must be satisfied. The formula is: in, This represents the number of points in the filtered point cloud. For point Vertical coordinates; This represents the actual terrain elevation. S1044: The optimized optical image data is used to screen the mapping area, which must simultaneously meet two quantitative indicators: sharpness and feature edge integrity. The clarity is measured by a structural similarity index. SSIM Assessment, requirements , Provide a clear reference image for the same area; the integrity of the ground feature edges is determined by edge matching degree. M Assessment, requirements ,in, , It is the overlap length between the extracted contour and the historical vector contour. It is the total length of the historical vector outline; Finally, the data of the surveyed area are integrated into... ,in, , For the first Optimized optical images For the first LiDAR point cloud after periodic filtering For the first The timestamp of the UAV remote sensing data collection period For the first Structural similarity index of optical images over time. For the first Edge matching degree of UAV remote sensing data in the period for Vertical error of terrain fitting during the period.

[0022] S2: Based on the aforementioned multi-period remote sensing dataset, associate the temporal characteristics of the same spatial location in the surveyed area, and through... NDVIThe vegetation growth rate is calculated by exponential time series fitting, and the land subsidence and subsidence trend are calculated by radar interferometry, outputting a dynamic change feature set.

[0023] For any spatial coordinate point within the surveying area Perform the following matching operation: S201: From the optimized optical image Extract coordinates Corresponding pixel grayscale value The focus is on screening the near-infrared band ( ) and red band ( The grayscale value of ) is then obtained from the first From the filtered lidar point cloud, obtain the coordinates ( Mean vertical coordinates of the surrounding effective point cloud , as coordinates No. The average topographic elevation over a period of time is calculated using the following formula: in, coordinates Number of effective point clouds in the surrounding area; For the first In point clouds, spatial coordinates fall on A collection of point clouds within the specified range; For point clouds The vertical coordinates; To obtain coordinates In the Time series datasets in the data : The process for calculating the vegetation growth rate is as follows: S202: Calculate the grayscale values ​​of the near-infrared and red bands of the optimized optical image. Period coordinates of NDVI The value, the formula is: in, For the first Period coordinates Near-infrared band ( The optimized grayscale value; For the first Period coordinates Red light band ( The optimized grayscale value; For the first Period coordinates The normalized vegetation index at the location, with a value range of 100%. Negative values ​​indicate non-vegetated areas (such as water bodies and bare rocks), while larger positive values ​​indicate more vigorous vegetation growth. Then, time difference conversion was performed, and the timestamps of each collection period were converted. (Formatted as "year-month-day") Convert to time difference based on period 1 (Unit: Quarter), the formula is: in, and Corresponding to the number in sequence The year and month in which the UAV remote sensing data was collected; and The data is sequentially assigned to the year and month of the first phase of UAV remote sensing data collection to ensure that the time difference is consistent with the time unit of the subsequent quarterly growth rate. With the time difference The independent variable is the normalized vegetation index. Using the dependent variable, construct a linear fitting model: in, For any time difference Below NDVI Predicted value; The slope of the fitted line, i.e., the coordinates Quarterly growth rate of vegetation at the location; The intercept of the fitted line represents the first period ( )of NDVI Initial value; Minimize the fitting error using the least squares method. Solve and The formula is: Validation is performed by calculating the goodness of fit. ,Require (indicating timing) NDVI The data fits the linear model well, and the growth rate calculation is reliable. The formula is: in, for Expect NDVI The mean, and ; For the first Expect NDVI Time series predicted values; The process for calculating the surface subsidence and subsidence trend is as follows: S203: Based on the topographic elevation of Phase 1 Based on this, calculate the first... Period coordinates Cumulative settlement (Negative values ​​represent settlement, positive values ​​represent uplift), the formula is: in, For the first Period coordinates The average topographic elevation; This is the baseline elevation value for Phase 1; time difference (Unit: Quarter) Further converted to time difference in "years" The formula is: This is to ensure that the time unit is consistent with the subsequent "annual settlement rate"; Subsequently For independent variable, Using the variable as the dependent variable, construct a linear fitting model and solve for the annual settlement rate. Negative values ​​represent average annual subsidence, and positive values ​​represent average annual uplift. The formula is: in, For any annual time difference, the cumulative settlement is predicted. The intercept of the fitted line represents the first period ( The initial settlement; Similarly, the solution is obtained using the least squares method. and The formula is: Finally, the effectiveness was verified, and the goodness of fit of the settlement was calculated. ,Require (Considering minor topographical fluctuations, the goodness-of-fit threshold is slightly lower than the vegetation growth rate, indicating a good fit between the time-series subsidence data and the linear model, and reliable calculation results for subsidence rate and trend). The formula is: in, for The average cumulative settlement over the period, and ; For the first Predicted cumulative settlement over a period of time; The integration result is the dynamically changing feature set, and the process is as follows: S204: For all valid coordinate points within the surveyed area It must meet the following conditions simultaneously. , as well as Construct a dynamically changing feature set ,and , Represents the effective mapping area; where the feature set of a single coordinate point is... Defined as: in, This refers to the quarterly growth rate of vegetation. The annual settlement rate; for NDVI Time series predicted values; This represents the predicted cumulative settlement value. This is the timestamp for the latest data collection.

[0024] S3: Based on the dynamic change feature set, preset the normal fluctuation range of vegetation growth rate and the safety threshold of land subsidence and verify them to determine whether it is necessary to trigger the manual interface. After manual confirmation, update the feature threshold or supplement the building feature library and output an automated surveying and analysis report.

[0025] For dynamically changing feature sets Each valid coordinate point The feature parameters are used to perform threshold checks on vegetation growth rate and land subsidence rate, respectively, to output the results of exceeding the limits. The process is as follows: S301: The vegetation growth rate threshold verification process is as follows: The normal fluctuation range of vegetation growth rate is set as follows: (Quarterly), recorded as ,in ; Actual quarterly vegetation growth rate at coordinate points The percentage difference from the normal range is used to determine whether the limit is exceeded. The formula is as follows: in, Relative deviation of vegetation growth rate (unitless); coordinate point The quarterly growth rate of vegetation; Standard quarterly growth rate for the corresponding vegetation type (unitless, based on a pre-defined vegetation type database, such as trees). Herbs ); like Marked as , indicates that the limit has been exceeded; if Marked as This indicates that it is normal. S302: The surface subsidence rate threshold verification process is as follows: The safety threshold for surface subsidence is set at 3. mm / The year was recorded as Negative values ​​indicate settlement, and the larger the absolute value, the more severe the settlement. Then, the absolute value of the actual annual settlement rate at the coordinate point is compared with the safety threshold to determine whether the limit has been exceeded. The formula is as follows: in, The surface subsidence rate exceeds the limit (unit: Year); coordinate point The annual settlement rate, with negative values ​​indicating settlement; Surface subsidence safety threshold (unit: Year); like Marked as , indicates that the limit has been exceeded; if Marked as This indicates that it is normal; S303: For each coordinate point Define the comprehensive overlimit flag as Flag The formula is: in, Represents coordinate points At least one parameter is out of limit. This indicates that both parameters are normal; The process for identifying unmatched man-made structures is as follows: S304: Based on the optimized optical image Extract coordinate points Surrounding building feature parameters Including the complexity of building outline Building height and spectral reflectance characteristics ; The complexity of the building outline The ratio of the number of contour inflection points to the contour length is calculated using an edge detection algorithm; the building height... The formula for calculating the vertical coordinate difference based on lidar point cloud is as follows: ,in coordinate point Maximum vertical coordinates of the point cloud in the surrounding area coordinate point Minimum vertical coordinates of the point cloud in the surrounding area (unit: The spectral reflectance characteristics Extracting common building materials (such as concrete and metal) in the near-infrared band ( The average reflectance of ) Subsequently, a comparison with the building feature database is performed, and a pre-set artificial building feature database is established. It includes the feature threshold range of common man-made structures, and ;in, The threshold range for building outline complexity. To minimize the contour complexity, a value of 0.05 is used. The maximum contour complexity is set to 0.3. For building height threshold range, The minimum building height is set at 2.5. m， The maximum building height is set to 50. m ; The minimum reflectivity is set to 0.3. The maximum reflectivity is set to 0.8. Define the matching degree as The formula is: in, Corresponding in sequence , and The Item building characteristic parameters; 、 For the first feature in the feature library Threshold range for each feature; This is an indicator function that returns 1 if the condition is met, and 0 otherwise. like That is, at least two features do not match, and are marked as , indicates no matching artificial structure; otherwise, it is marked as , indicating that a match has been found; When the coordinates satisfy or When the manual intervention interface is triggered, the relevant parameters are updated after calibration, and the final verification result is output. The process is as follows: S305: The intervention trigger condition is all valid coordinate points within the survey area. or The point triggers a manual intervention interface with feature annotation, and the annotation content includes the point exceeding the limit. , Calculation results and unmatched buildings Details of feature parameter comparison with feature library; Subsequently, parameter calibration and updates were performed, including vegetation growth rate threshold calibration, land subsidence threshold calibration, and building feature library supplementation: The process of calibrating the vegetation growth rate threshold involves manually adjusting the standard growth rate according to the actual vegetation type. The updated threshold range is denoted as The formula is , ,in (Fixed fluctuation ratio) The calibrated standard growth rate; The process for calibrating the surface subsidence threshold is as follows: if it is manually confirmed that a higher subsidence rate is permissible in a local area, the safety threshold is updated. And record the calibration instructions; The process of supplementing the building feature database involves adding unmatched areas that are identified as man-made buildings. Parameters added to the building feature library Update threshold range ; S306: Based on the calibrated parameters, recalculate the comprehensive verification result for each coordinate point. The formula is: in, This indicates that it is normal; This indicates that the vegetation growth rate has exceeded the limit; This indicates that the surface subsidence has exceeded the limit; This indicates that both parameters have exceeded the limit; This indicates that no matching artificial structure was found. 、 , and All are recalculated marker values ​​after calibration; S307: Integrate all the above calculation results and verification information to generate an automated surveying and mapping analysis report. The automated surveying and mapping analysis report includes an overview of the surveying area, dynamic change statistics, verification result details, and recommended measures. The survey area overview includes the effective survey area range, data collection period, and data quality indicators for each period (clarity, feature edge integrity); the dynamic change statistics include the mean and standard deviation of vegetation growth rate, the percentage of areas exceeding limits, the mean and standard deviation of surface subsidence rate, the percentage of areas exceeding limits, and the number and distribution of mismatched artificial structures; the verification result details include each coordinate point. The key parameter calculation values ​​and calibration records; the recommended measures propose monitoring frequency adjustment suggestions for areas exceeding the limits, and feature library update suggestions for buildings that do not match.

[0026] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0027] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. An automated mapping and analysis method based on UAV remote sensing data, characterized in that, include: S1: Acquire multiple periods of UAV remote sensing data and perform scenario-based UAV remote sensing data optimization processing, extract the mapping areas that meet the conditions in each period of data, and output the multi-period remote sensing dataset; S2: Based on the aforementioned multi-period remote sensing dataset, associate the temporal characteristics of the same spatial location in the surveyed area, and through... NDVI Exponential time series fitting is used to calculate vegetation growth rate, radar interferometry is used to calculate land subsidence and subsidence trend, and a dynamic change feature set is output. S3: Based on the dynamic change feature set, preset the normal fluctuation range of vegetation growth rate and the safety threshold of land subsidence and verify them to determine whether it is necessary to trigger the manual interface. After manual confirmation, update the feature threshold or supplement the building feature library and output an automated surveying and analysis report.

2. The automated mapping and analysis method based on UAV remote sensing data according to claim 1, characterized in that, The following data was acquired from multiple UAV remote sensing datasets: S101: The acquired multi-period UAV remote sensing data set is... ;in, This represents the total number of data collection periods, and ; For the first The first phase of the UAV remote sensing dataset includes the first phase. A subset of optical images from the period LiDAR point cloud subset and UAV remote sensing data acquisition timestamp .

3. The automated mapping and analysis method based on UAV remote sensing data according to claim 1, characterized in that, The workflow for optimizing and processing drone remote sensing data for different scenarios is as follows: S102: For optical images with strong backlighting, an optimized optical image is obtained using a multi-exposure fusion algorithm, with the following formula: in, These are the planar coordinates of the pixel. This represents the total number of optical image frames. For optical image frame indexing; To optimize the optical image at the pixel level The grayscale value at that location; For the first Issue No. Frame optical images in The original grayscale value of a pixel; For the first Frame optical images in The weights of the pixels satisfy the conditions ; S103: To correct the reflectivity deviation in flooded areas after heavy rain, the dark target method is used. The formula is: in, To correct the water body at wavelength Reflectance at that location; To correct the water body at wavelength The original reflectivity at that location is determined by of DN Value conversion, and , wavelength Image at the location DN value, wavelength The saturated radiance at that location For corresponding saturation DN value; For dark targets at wavelength Reflectance at that location; S104: For isolated noise points in the lidar point cloud, the processing steps are as follows: S1041: Calculate the subset of the lidar point cloud The Middle Points The set of neighborhood points is given by the formula: in, for Surrounding distance ≤ The set of all neighboring points; The neighborhood search radius; for The Middle One point; For point 3D coordinates; S1042: Based on the neighborhood point set Calculation points Average distance to all points in the neighborhood To reflect the point The degree of fit with the surrounding terrain was determined, and the average distance of all points in the neighborhood was calculated. mean and standard deviation ; S1043: Based on the average distance of all points in the neighborhood mean and standard deviation Noise point identification and removal: If Judgment point These are noise points, and they should be removed; if Judgment point Points that are considered valid are retained; the final retained valid points constitute a point cloud. The terrain fitting vertical error must be satisfied. The formula is: in, This represents the number of points in the filtered point cloud. For point Vertical coordinates; This represents the actual terrain elevation. S1044: The optimized optical image data is used to screen the mapping area, which must simultaneously meet two quantitative indicators: sharpness and feature edge integrity. The clarity is measured by a structural similarity index. SSIM Assessment, requirements , Provide a clear reference image for the same area; the integrity of the ground feature edges is determined by edge matching degree. M Assessment, requirements ; Finally, the data of the surveyed area are integrated into... ,in, , corresponding to the number The optimized optical image, filtered lidar point cloud, UAV remote sensing data acquisition timestamp, optical image structure similarity index, UAV remote sensing data edge matching degree, and terrain fitting vertical error are analyzed.

4. The automated mapping and analysis method based on UAV remote sensing data according to claim 3, characterized in that, For any spatial coordinate point within the surveying area Perform the following matching operation: S201: From the optimized optical image Extract coordinates Corresponding pixel grayscale value The focus is on filtering the grayscale values ​​in the near-infrared and red light bands, and then obtaining the coordinates from the filtered lidar point cloud. Mean vertical coordinates of the surrounding effective point cloud , as coordinates No. The average topographic elevation over a period of time; to obtain coordinates. In the Time series datasets in the data : The process for calculating the vegetation growth rate is as follows: S202: Calculate the grayscale values ​​of the near-infrared and red bands of the optimized optical image. Period coordinates of NDVI The value, the formula is: in, For the first Period coordinates Optimized grayscale values ​​in the near-infrared band; For the first Period coordinates Optimized grayscale values ​​in the red band; For the first Period coordinates The normalized vegetation index at the location, with a value range of Negative values ​​indicate non-vegetated areas, while larger positive values ​​indicate more vigorous vegetation growth. Then, time difference conversion was performed, and the collection timestamps for each period were formatted as year, month, and day. Converted to time difference based on period 1 The unit is quarters, based on the time difference mentioned above. The independent variable is the normalized vegetation index. Using the dependent variable, construct a linear fitting model: in, For any time difference Below NDVI Predicted value; The slope of the fitted line, i.e., the coordinates Quarterly growth rate of vegetation at the location; The intercept of the fitted line represents the first period. of NDVI Initial value; Minimize the fitting error using the least squares method. Solve and And by calculating the goodness of fit To perform validity verification, the following requirements must be met. to indicate time sequence NDVI The data fits the linear model well, and the growth rate calculation results are reliable.

5. The automated mapping and analysis method based on UAV remote sensing data according to claim 4, characterized in that, The process for calculating the surface subsidence and subsidence trend is as follows: S203: Based on the topographic elevation of Phase 1 Based on this, calculate the first... Period coordinates Cumulative settlement Negative values ​​represent settlement, and positive values ​​represent uplift. The formula is: in, For the first Period coordinates The average topographic elevation; This is the baseline elevation value for Phase 1; Time difference in quarters Further converted into time difference in years ; Subsequently For independent variable, Using the variable as the dependent variable, construct a linear fitting model and solve for the annual settlement rate. Negative values ​​represent average annual subsidence, and positive values ​​represent average annual uplift. The formula is: in, For any annual time difference, the cumulative settlement is predicted. The intercept of the fitted line represents the first period. The initial settlement; Similarly, the solution is obtained using the least squares method. and And by calculating the goodness of fit of the settlement To perform validity verification, the following requirements must be met. This indicates that the time-series settlement data fits the linear model well, and the calculation results of settlement rate and trend are reliable; The integration result is the dynamically changing feature set, and the process is as follows: S204: For all valid coordinate points within the surveyed area It must meet the following conditions simultaneously. , as well as Construct a dynamically changing feature set ,and , Represents the effective mapping area; where the feature set of a single coordinate point is... Defined as , This is the timestamp for the latest data collection.

6. The automated mapping and analysis method based on UAV remote sensing data according to claim 5, characterized in that, For dynamically changing feature sets Each valid coordinate point The feature parameters are used to perform threshold checks on vegetation growth rate and land subsidence rate, respectively, to output the results of exceeding the limits. The process is as follows: S301: The vegetation growth rate threshold verification process is as follows: The vegetation growth rate threshold verification is performed using the actual quarterly vegetation growth rate at coordinate points. Relative deviation from the normal range Implementation; if Marked as , indicates that the limit has been exceeded; if Marked as , indicating normal; among them and All are labeled variables; S302: The surface subsidence rate threshold verification process is as follows: The surface subsidence rate threshold verification is performed by checking the coordinate points. The absolute value of the actual annual settlement rate is compared with the safety threshold for exceeding the limit. Calculation implementation; if Marked as , indicates that the limit has been exceeded; if Marked as This indicates normal; similarly, it includes... To mark variables; S303: For each coordinate point Define the comprehensive overlimit flag as Flag The formula is: in, Represents coordinate points At least one parameter is out of limit. This indicates that both parameters are normal; After completing the threshold verification of vegetation growth rate and land subsidence rate, the identification of unmatched artificial structures is carried out, and the process is as follows: S304: Based on the optimized optical image Extract coordinate points Surrounding building feature parameters Including the complexity of building outline Building height and spectral reflectance characteristics ; Subsequently, a comparison with the building feature database is performed, and a pre-set artificial building feature database is established. It includes the feature threshold range of common man-made structures, and ;in, This refers to the threshold range for building outline complexity. To minimize the contour complexity, The maximum contour complexity; For building height threshold range, Minimum building height ， Maximum building height; Minimum reflectivity Maximum reflectivity; Define the matching degree as ,like That is, at least two features do not match, and are marked as , indicates no matching artificial structure; otherwise, it is marked as , indicates that a match has been found; similarly, in Used to label variables.

7. The automated mapping and analysis method based on UAV remote sensing data according to claim 6, characterized in that, When coordinate point satisfy or When the manual intervention interface is triggered, the relevant parameters are updated after calibration, and the final verification result is output. The process is as follows: S305: The intervention trigger condition is all valid coordinate points within the survey area. or The point triggers a manual intervention interface with feature annotation, and the annotation content includes the point exceeding the limit. , Calculation results and unmatched buildings Details of feature parameter comparison with feature library; Subsequently, parameter calibration and updates were performed, including vegetation growth rate threshold calibration, land subsidence threshold calibration, and building feature library supplementation. S306: Based on the calibrated parameters, recalculate the comprehensive verification result for each coordinate point. The formula is: in, This indicates that it is normal; This indicates that the vegetation growth rate has exceeded the limit; This indicates that the surface subsidence has exceeded the limit; This indicates that both parameters have exceeded the limit; This indicates that no matching artificial structure was found. 、 , and All are recalculated marker values ​​after calibration; S307: Integrate all calculation results and verification information to generate an automated surveying and mapping analysis report. The automated surveying and mapping analysis report includes an overview of the surveying area, dynamic change statistics, verification result details, and recommended measures.

8. An automated mapping and analysis system based on UAV remote sensing data, wherein the system is applied to the automated mapping and analysis method based on UAV remote sensing data as described in claims 1-7, characterized in that, The system includes the following steps: M1: Data acquisition module, which acquires remote sensing data from multiple periods through the UAV data acquisition unit and transmits it to the scene-specific optimization processing unit; it performs scene-specific processing such as multi-exposure fusion optimization, dark target reflectivity correction, and neighborhood analysis noise removal for strong backlight optical images, post-rainstorm water accumulation area images, and lidar point clouds respectively; then, through the region filtering unit, it extracts the mapping areas in each period of data that meet the requirements of sharpness ≥ 0.9 and ground feature edge integrity ≥ 0.85, and outputs multi-period remote sensing datasets; M2: Dynamic feature calculation module, which uses the time-series matching unit to associate time-series information of the same spatial location in the surveyed area from multiple remote sensing datasets to construct a coordinate point time-series dataset; the vegetation rate calculation unit is based on... NDVI The exponential time series fitting model solves and outputs the quarterly vegetation growth rate; the subsidence trend calculation unit uses radar interferometry technology to calculate the cumulative surface subsidence and annual subsidence rate, integrates growth rate and subsidence trend data to form a dynamic change feature set and outputs it. M3: Feature verification and report generation module. The threshold verification unit presets the normal fluctuation range of vegetation growth rate and the safety threshold of land subsidence, and judges the parameters in the dynamic feature set that exceed the limit. The building identification unit extracts the building features around the coordinate point and compares them with the preset building feature library to identify unmatched artificial buildings. If there are buildings with parameters exceeding the limit or unmatched buildings, the manual intervention interface is triggered. After manual confirmation, the feature threshold is updated or the building feature library is supplemented. Finally, the report integration unit outputs an automated surveying and analysis report.

Images