Unmanned aerial vehicle remote sensing Gaussian splashing change detection method for dynamic monitoring of natural resources

Abstract

The invention discloses a natural resource dynamic monitoring-oriented unmanned aerial vehicle remote sensing Gaussian splashing change detection method, particularly relates to the field of natural resource monitoring, and aims to solve the problem that natural evolution and man-made interference are difficult to distinguish. Separating natural and artificial change characteristics by combining time axis gradient analysis, predicting a density field in a non-interference state by using a low-frequency gradient component, accurately positioning an artificial interference area through multi-scale difference fusion and dynamic threshold setting, and finally constructing a causal reasoning network to reveal an interference source and a propagation path; meanwhile, by fusing geographic space data and historical activity information, causal association analysis of interference events is realized, a targeted basis is provided for ecological protection decision, timeliness of natural resource supervision is effectively improved, and efficient support is provided for resource management and environmental protection.

Description

Technical Field

[0001] The present invention relates to the field of natural resource monitoring, and more specifically, to a method for detecting Gaussian splash changes using unmanned aerial vehicle remote sensing for dynamic monitoring of natural resources. Background Art

[0002] In the field of dynamic monitoring of natural resources, existing technologies generally adopt density difference detection methods based on fixed thresholds. The core limitation of existing technologies is that they cannot effectively model the essential differences between natural evolution processes and human intervention behaviors. Natural evolution usually manifests itself as gradual changes that conform to ecological laws over a long time scale, while human intervention often causes local mutations that violate natural laws in a short period of time. The two may show similar characteristics in the amplitude of density field changes, but the dynamic patterns in the time dimension are completely different. Due to the lack of quantitative modeling capabilities for the rate of temporal changes, existing methods mistakenly classify slow natural processes and sudden human events as similar changes; at the same time, external environmental noise and sensor acquisition errors form pseudo-change signals in the density field, further exacerbating the confusion of detection results. This defect makes it difficult for the monitoring system to accurately distinguish key scenes such as seasonal vegetation withering and flourishing and natural river diversion, resulting in the dual problems of increased false alarm rate and missed detection of key events. More seriously, existing technologies cannot parse the causal relationship between change events. For example, it is impossible to identify whether vegetation degradation in a certain area is indirectly caused by neighboring mining activities, resulting in a lack of pertinence in subsequent governance measures. In order to solve the above problems, a technical solution is now provided. Summary of the invention

[0003] In order to overcome the above-mentioned defects of the prior art, an embodiment of the present invention provides a UAV remote sensing Gaussian splash change detection method for dynamic monitoring of natural resources. It generates a spatiotemporal continuous observation density field through multi-phase point cloud data conversion, separates natural and man-made change characteristics in combination with time axis gradient analysis, and uses low-frequency gradient components to predict the density field under interference-free conditions. It also accurately locates the man-made interference area through multi-scale difference fusion and dynamic threshold setting, and finally constructs a causal reasoning network to reveal the source of interference and the propagation path. At the same time, by integrating geographic spatial data with historical activity information, it realizes the causal correlation analysis of interference events, provides targeted basis for ecological protection decision-making, and provides efficient support for resource management and environmental protection, so as to solve the problems raised in the above-mentioned background technology.

[0004] To achieve the above object, the present invention provides the following technical solutions:

[0005] S1. Obtain multi-temporal UAV remote sensing point cloud data and convert it into a spatiotemporal continuous observation density field;

[0006] S2. Perform time axis gradient analysis on the observed density field, separate the low frequency gradient components of natural evolution and the high frequency gradient components of artificial mutation, and generate a classification map of change patterns;

[0007] S3. Based on the naturally evolving low-frequency gradient component, obtain the predicted density field under the undisturbed state;

[0008] S4. Calculate the difference between the observed density field and the predicted density field at different spatial scales, and fuse the differences between scales into comprehensive differences through consistency weighted fusion. Combine the historical change amplitude of the region and the spatial heterogeneity to dynamically set the judgment threshold, and locate the human interference area based on the comparison between the comprehensive difference and the dynamic threshold.

[0009] S5. Combine geospatial data and historical activity information to construct a causal reasoning network to analyze the causal sources and propagation paths of human disturbance areas.

[0010] In a preferred embodiment, step S1 includes the following contents:

[0011] The multi-temporal UAV remote sensing point cloud data are denoised, and a density-based clustering method is used to remove outliers to generate denoised point cloud data. All denoised point cloud data are aligned to a unified geographic coordinate system through a registration method. The study area is divided into a uniform two-dimensional grid and subdivided into multiple altitude layers in the vertical direction. The denoised point cloud data are assigned to the corresponding two-dimensional grid cells and altitude layers according to the three-dimensional coordinates. The local observation density value is calculated for the point set in each altitude layer using the Gaussian kernel density estimation method, and the comprehensive observation density value is generated by weighted summation. The comprehensive observation density value is linearly interpolated in the time dimension to fill the missing phase data, and Gaussian filtering is applied for spatiotemporal smoothing to finally generate the observation density field.

[0012] In a preferred embodiment, step S2 includes the following contents:

[0013] The observed density field data is subjected to time axis gradient analysis. The density change gradient is calculated by the central difference method for the observed density value sequence of each spatial grid unit. The density change gradient sequence is converted into frequency domain spectrum data by discrete Fourier transform. The spectrum data is separated into low-frequency components and high-frequency components according to a preset frequency threshold. The low-frequency components and high-frequency components are respectively reconstructed into low-frequency gradient components and high-frequency gradient components in the time domain by inverse discrete Fourier transform. The relative intensity ratio of the high-frequency gradient component to the low-frequency gradient component is calculated and compared with the classification threshold to generate a change pattern classification map. Finally, the low-frequency gradient component field, high-frequency gradient component field and change pattern classification map are output as input for predicting the density field and locating the human interference area.

[0014] In a preferred embodiment, step S3 includes the following contents:

[0015] Based on the separated low-frequency gradient components, the observed density values ​​at the initial time points provided by the observed density field are used to reconstruct the predicted density values ​​at subsequent time points by recursively accumulating the product of the low-frequency gradient components and the time interval. The reconstructed predicted density field is smoothed by applying Gaussian filtering in the time dimension. The smoothed predicted density values ​​are generated by taking Gaussian weighted average of the predicted density values ​​at adjacent time points in the time window. The smoothed predicted density values ​​are integrated into a time-space three-dimensional tensor to form a predicted density field in an undisturbed state.

[0016] In a preferred embodiment, step S4 includes the following contents:

[0017] For the observed density field and the predicted density field, the absolute value of the difference between the smoothed observed density value and the smoothed predicted density value is calculated by Gaussian smoothing at multiple spatial scales as the difference value. The consistent weighted fusion method is used to calculate the entropy value of the difference field in the local neighborhood at each spatial scale, and the weight is determined according to the entropy value. The comprehensive difference value is generated by weighted average. The dynamic threshold is set based on the spatial heterogeneity calculated by combining the average absolute difference of the historical observed density value changes of each spatial grid unit and the local Moran's I index. The map of the human interference area is generated by comparing the comprehensive difference value with the dynamic threshold.

[0018] In a preferred embodiment, step S5 includes the following contents:

[0019] The human interference area map, geospatial data and historical activity information are integrated and aligned in the time and space dimensions. After defining the interference source nodes based on the historical activity information and the interference area nodes based on the human interference area map, a causal network containing temporal causal edges, spatial causal edges and source-area edges is constructed. Based on the screened causal network, the shortest path algorithm is used to identify the interference sources and their propagation paths, and the main interference sources and the impact range are counted.

[0020] In a preferred embodiment, step S5 further includes the following contents:

[0021] The temporal correlation of nodes in the interference region at adjacent time points is calculated using the Pearson correlation coefficient and compared with the time threshold to screen the temporal causal edges.

[0022] In a preferred embodiment, step S5 further includes the following contents:

[0023] The spatial correlation between adjacent grids at the same time point is calculated by the Moran's I index and compared with the spatial threshold to screen the spatial causal edges.

[0024] In a preferred embodiment, step S5 further includes the following contents:

[0025] The causal impact of the interference source node on the interference region node is calculated through the Granger causality test and compared with the significance level to screen the source-region edge.

[0026] The technical effects and advantages of the unmanned aerial vehicle remote sensing Gaussian splash change detection method for dynamic monitoring of natural resources of the present invention are as follows:

[0027] The present invention significantly improves the accuracy and scientificity of natural resource change monitoring through the unmanned aerial vehicle remote sensing Gaussian splash change detection method for dynamic monitoring of natural resources. In view of the limitations of the fixed threshold method in the prior art that it is difficult to distinguish between natural evolution and human intervention, is susceptible to noise interference, and lacks causal analysis, the present invention generates a spatiotemporal continuous observation density field through multi-phase point cloud data conversion, combines time axis gradient analysis to separate natural and human change characteristics, uses low-frequency gradient components to predict the density field under interference-free conditions, and accurately locates the human interference area through multi-scale difference fusion and dynamic threshold setting, and finally constructs a causal reasoning network to reveal the source of interference and propagation path. This technical solution overcomes the defects of traditional methods in insufficient modeling of time series change rate and confusion of pseudo-change signals, reduces false alarm rate and missed detection rate, and enhances the ability to distinguish complex scenes. At the same time, by integrating geospatial data with historical activity information, the present invention realizes the causal correlation analysis of interference events, providing a targeted basis for ecological protection decision-making. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 This is a flow chart of the Gaussian splash change detection method for unmanned aerial vehicle remote sensing for dynamic monitoring of natural resources according to the present invention. DETAILED DESCRIPTION

[0029] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0030] Embodiment 1: Figure 1 The present invention provides a method for detecting Gaussian splash changes using remote sensing by an unmanned aerial vehicle for dynamic monitoring of natural resources, comprising:

[0031] S1. Obtain multi-temporal UAV remote sensing point cloud data and convert it into a spatiotemporally continuous observation density field.

[0032] S2. Perform time axis gradient analysis on the observed density field to separate the low-frequency gradient components of natural evolution and the high-frequency gradient components of artificial mutations, and generate a classification map of change patterns.

[0033] S3. Based on the naturally evolving low-frequency gradient component, the predicted density field under undisturbed conditions is obtained.

[0034] S4. Calculate the difference between the observed density field and the predicted density field at different spatial scales, and fuse them into comprehensive differences through weighted fusion of the consistency of the differences between scales. Dynamically set the judgment threshold based on the historical change amplitude of the region and spatial heterogeneity, and locate the human interference area based on the comparison between the comprehensive difference and the dynamic threshold.

[0035] S5. Combine geospatial data and historical activity information to construct a causal reasoning network to analyze the causal sources and propagation paths of human disturbance areas.

[0036] In the field of dynamic monitoring of natural resources, real-time grasp of phenomena such as changes in vegetation cover, river morphology adjustments, and changes in land use patterns is crucial for ecological protection and resource management. These changes are often jointly affected by natural evolution processes and human intervention behaviors. To this end, it is necessary to obtain high-resolution multi-phase point cloud data through UAV remote sensing technology to construct observation density field data that can reflect the spatiotemporal dynamics of surface characteristics, thereby providing a data basis for the subsequent distinction between natural and human changes, locating interference areas, and analyzing causal relationships. Step S1 generates spatiotemporally continuous observation density field data by processing multi-phase UAV remote sensing point cloud data, laying the foundation for the implementation of the entire monitoring method.

[0037] Step S1 includes the following contents:

[0038] S1.1, Data acquisition:

[0039] First, obtain multi-temporal UAV remote sensing point cloud data. These data are point cloud collections collected at different time points. Each set of point cloud data contains a large number of three-dimensional coordinate information of points, namely the horizontal coordinate, the horizontal coordinate, and the vertical coordinate. These multi-temporal UAV remote sensing point cloud data are the basic input data for subsequent observation density field conversion.

[0040] Multi-temporal UAV remote sensing point cloud data refers to a collection of point cloud data collected at multiple different time points by remote sensing equipment carried by UAVs. Specifically, remote sensing equipment (such as lidar or high-resolution cameras) repeatedly acquires three-dimensional spatial information of the surface at different time periods (such as daily, monthly or annually) for the same study area during the flight of the UAV. This information is recorded in the form of a point cloud, and each point contains a horizontal coordinate, a horizontal coordinate, and a vertical coordinate, and may also include additional attributes such as reflection intensity or color value. The so-called multi-temporal emphasizes that data collection has time series characteristics and can reflect the dynamic process of surface characteristics changing over time, such as vegetation growth, river diversion or human development.

[0041] S1.2, point cloud preprocessing:

[0042] Each set of multi-temporal UAV remote sensing point cloud data is denoised by using a density-based clustering method. By analyzing the distance between points, outliers that deviate from the main point cloud distribution are identified and removed to generate denoised point cloud data. After denoising, the registration method is used to align all denoised point cloud data into a unified geographic coordinate system by comparing the spatial position differences between point cloud data collected at different time points, ensuring the consistency of the multi-temporal point cloud data in spatial position.

[0043] S1.3, spatial grid division:

[0044] The study area is divided into a uniform two-dimensional grid, and the horizontal and vertical spacing of each grid unit is a pre-set fixed value, which is used to represent a specific spatial position. On this basis, each grid unit is further subdivided into multiple height layers along the vertical direction, that is, the elevation direction, and the interval between each height layer is a pre-set fixed value, forming a three-dimensional space division structure.

[0045] Two-dimensional grid division discretizes the continuous space of the study area into calculable units, which is convenient for subsequent density calculation and spatial distribution analysis. Height layer subdivision can reflect the distribution differences of point cloud data in the vertical direction, such as distinguishing vegetation canopy from surface features, and enhance the comprehensiveness of spatial description.

[0046] S1.4, Observed density field calculation:

[0047] For each set of denoised multi-temporal UAV remote sensing point cloud data, all points are assigned to the corresponding two-dimensional grid cells and their altitude layers according to the three-dimensional coordinates of the points. For each grid cell point set in a specific altitude layer, the Gaussian kernel density estimation method is used to calculate the local density. The specific calculation process is as follows: first determine the distance between each point in the altitude layer and the geometric center of the grid cell, then assign a weight to each point according to the distance through the Gaussian function, where the weight size is controlled by the preset smoothing bandwidth parameter; then, add the weight values ​​of all points in the altitude layer and perform normalization to obtain the local observation density value of the altitude layer. After that, the local observation density values ​​of each grid cell in all altitude layers are weighted and summed to generate the comprehensive observation density value of the grid cell at the current time point, where the weight coefficient is pre-adjusted according to the type of ground object (such as vegetation or building).

[0048] The Gaussian kernel density estimation method reduces the local noise in the point cloud data distribution through smoothing and provides a continuous and stable density representation. The weighted integration of the local density of the altitude layer comprehensively considers the distribution characteristics of the ground objects in the vertical direction and enhances the ability of the comprehensive observation density value to characterize the surface characteristics.

[0049] S1.5, spatiotemporal continuity processing:

[0050] For missing time points in multi-temporal UAV remote sensing point cloud data, the comprehensive observation density values ​​of the adjacent time points are used to calculate the comprehensive observation density values ​​of the missing time points through linear interpolation to fill the gaps in the time series. After the interpolation is completed, the Gaussian filtering method is applied to the comprehensive observation density values ​​of all time points and spatial positions, and smoothing is performed in the time dimension and space dimension to generate smoothed observation density field data.

[0051] The linear interpolation method infers missing values ​​through the trend of adjacent data, ensuring the integrity of the time series and supporting subsequent analysis based on the time dimension. The smoothing process of Gaussian filtering in the time and space dimensions reduces the impact of short-term fluctuations and spatial noise, and improves the stability and continuity of the observed density field data.

[0052] S1.6, Observed density field output:

[0053] Finally, a spatiotemporally continuous observation density field data structure is generated, which contains the smoothed integrated observation density values ​​of all time points and all grid cells, forming a three-dimensional data set whose dimensions are time, number of grid rows, and number of grid columns for subsequent analysis.

[0054] Step S1 has generated spatiotemporally continuous observation density field data by processing multi-temporal UAV remote sensing point cloud data, which fully reflects the distribution characteristics of surface features in time and space. Step S2 inherits this observation density field data and aims to separate the low-frequency gradient components of natural evolution and the high-frequency gradient components of artificial mutations through time axis gradient analysis, and generate a change pattern classification map to reveal the inherent laws and abnormal characteristics of natural resource changes.

[0055] Step S2 includes the following contents:

[0056] S2.1, time axis gradient calculation:

[0057] The observed density field data generated in step S1 is processed, and the change rate of the observed density value sequence of each spatial grid unit at different time points is calculated along the time dimension. The calculation method adopts the central difference method, that is, the difference between the observed density value at the previous time point and the observed density value at the next time point is divided by twice the time interval between the two time points to obtain the density change gradient of the corresponding spatial grid unit at the current time point. This density change gradient reflects the change speed of the observed density value in the time dimension, providing basic data for subsequent analysis.

[0058] The density gradient can directly quantify the trend of the observed density value over time. The natural evolution of natural resources usually manifests as slow and continuous changes, while mutations caused by human intervention present as rapid and drastic changes. By calculating the density gradient, these trends can be converted into quantifiable values, providing data support for frequency domain analysis and component separation.

[0059] S2.2, frequency domain conversion:

[0060] The density change gradient sequence of each spatial grid unit is converted into the frequency domain, and the discrete Fourier transform method is used to convert the density change gradient data in the time dimension into spectrum data in the frequency dimension. This conversion process decomposes the density change gradient sequence into multiple sinusoidal wave components of different frequencies, each sinusoidal wave component corresponds to a specific frequency value, reflecting the periodic characteristics of the density change gradient in the time dimension.

[0061] Frequency domain conversion can decompose the density change gradient in the time dimension into components of different frequencies, making it easier to distinguish the periodic characteristics of the change. Natural evolution usually manifests as slow changes in low frequencies, while artificial mutations usually manifest as rapid changes in high frequencies. Through this conversion, a clear frequency domain basis can be provided for the separation of low-frequency components and high-frequency components, improving the accuracy of the analysis.

[0062] S2.3, separation of low-frequency and high-frequency components:

[0063] In the frequency domain, the spectrum data is divided into low-frequency components and high-frequency components according to the pre-set frequency threshold. The determination of the frequency threshold is based on the time scale of natural resource evolution, and the maximum frequency value of natural changes is statistically analyzed by analyzing historical data. Spectral components below this frequency threshold are classified as low-frequency components, representing the characteristics of natural evolution; spectral components above this frequency threshold are classified as high-frequency components, representing the characteristics of artificial mutations.

[0064] By separating the spectral components through frequency thresholds, we can ensure that the low-frequency components accurately reflect the slow-changing trend of natural evolution, while the high-frequency components can effectively capture the rapid changes caused by human intervention. This separation method improves the pertinence of the analysis, allowing the characteristics of natural evolution and artificial mutation to be clearly distinguished in the frequency domain.

[0065] S2.4, inverse transformation and component reconstruction:

[0066] The separated low-frequency components and high-frequency components are respectively subjected to inverse discrete Fourier transform, and are converted from the frequency domain back to the time domain to generate the time series of low-frequency gradient components and the time series of high-frequency gradient components. This reconstruction process restores the change characteristics of low-frequency components and high-frequency components in the time dimension, which facilitates subsequent analysis and classification in the context of the time series.

[0067] The inverse discrete Fourier transform converts the low-frequency and high-frequency components in the frequency domain back to the time domain, so that these components can be intuitively understood and analyzed within the time frame of the original density change gradient sequence. This reconstruction method provides a data form that is easy to process for the subsequent classification of change patterns, ensuring the continuity and consistency of the analysis process.

[0068] S2.5, Classification of change patterns:

[0069] For each spatial grid unit at each time point, the relative intensity ratio of the high-frequency gradient component to the low-frequency gradient component is calculated. The specific calculation method is to take the ratio of the absolute value of the high-frequency gradient component to the absolute value of the low-frequency gradient component, and compare this ratio with the pre-set classification threshold. If this relative intensity ratio is greater than the classification threshold, it is determined that the change pattern of the spatial grid unit at this time point is an artificial mutation; if this relative intensity ratio is less than or equal to the classification threshold, it is determined to be a natural evolution.

[0070] The relative intensity ratio can quantify the degree of dominance of the high-frequency gradient component in the density change gradient, reflecting the contribution of artificial mutations to the total change.

[0071] S2.6, output results:

[0072] Finally, low-frequency gradient component field, high-frequency gradient component field and change pattern classification map are generated. The low-frequency gradient component field represents the natural evolution trend of each spatial grid unit in the time dimension. The high-frequency gradient component field and change pattern classification map are used to identify and locate the human intervention area, reflecting the distribution and intensity of human mutations.

[0073] Step S2 receives the observed density field data generated in step S1, and separates the low-frequency gradient components representing natural evolution and the high-frequency gradient components representing artificial mutations through a series of processes such as time axis gradient calculation, frequency domain conversion, separation of low-frequency and high-frequency components, inverse transformation and component reconstruction, and change pattern classification, and generates a change pattern classification map. These processing results provide key input data for step S3 to predict the observed density field under the undisturbed state and step S4 to locate the artificial intervention area, ensuring the logic and integrity of the entire method.

[0074] Step S2 successfully separates the naturally evolving low-frequency gradient components and the artificially mutated high-frequency gradient components by performing time-axis gradient calculation and frequency-domain processing on the observed density field data, and generates a change pattern classification map based on the relative intensity ratio.

[0075] Step S2 has separated the low-frequency gradient components of natural evolution and the high-frequency gradient components of artificial mutation by performing time axis gradient analysis on the spatiotemporal continuous observed density field data, and generated a classification map of change patterns, in which the low-frequency gradient components reflect the trend of natural changes in the observed density field over time and are not affected by human interference. Step S3 uses this low-frequency gradient component as input, and generates a predicted density field under interference-free conditions through reconstruction and smoothing, in order to provide benchmark data for step S4, so as to perform difference analysis with the observed density field and locate the area of ​​human interference.

[0076] Step S3 includes the following contents:

[0077] S3.1, using low-frequency gradient components to reconstruct the natural evolution trend:

[0078] In the process of generating the predicted density field under the interference-free state, it is first necessary to use the low-frequency gradient component to reconstruct the natural evolution trend of the observed density field in the time dimension. The low-frequency gradient component is obtained through the time axis gradient analysis in step S2, which represents the rate of natural change of the observed density field over time, and the influence of artificial mutations has been eliminated. The processing method is to start from the observed density value at the initial time point, gradually accumulate the low-frequency gradient components, and predict the observed density values ​​at subsequent time points. The specific operation steps are: starting from the observed density value at the initial time point provided by the observed density field generated in step S1, for each spatial grid unit, calculate the predicted density value at the subsequent time point. The predicted density value is calculated by adding the predicted density value at the previous time point to the product of the low-frequency gradient component and the time interval, and recursively obtaining the predicted density values ​​at all time points.

[0079] The low-frequency gradient component only reflects the trend of natural evolution and does not contain the mutation component of human interference. Therefore, the observed density field reconstructed based on the low-frequency gradient component can truly simulate the natural changes under the undisturbed state. The recursive accumulation calculation method is simple and efficient, and can continuously predict the density value in the time series to ensure the consistency of the prediction results with the natural evolution trend.

[0080] S3.2, smoothing for improved stability:

[0081] After reconstructing the natural evolution trend using low-frequency gradient components, the directly obtained predicted density field needs to be further smoothed due to the possible presence of tiny noise or local discontinuity in the low-frequency gradient components to improve its spatiotemporal continuity and reliability. The processing method is to apply Gaussian filtering to the reconstructed predicted density field in the time dimension, and generate a smoothed predicted density field by weighted averaging the predicted density values ​​of adjacent time points. Specifically, a fixed time window is selected, and the predicted density values ​​of each time point are weighted and summed within this window. The weighting method uses the weight of the Gaussian distribution. The weight value is determined according to the distance between the time point and the center point. The closer the distance, the higher the weight. The sum of all weights is 1, and finally a smooth predicted density value is obtained.

[0082] S3.3, generate predicted density field:

[0083] After the smoothing process is completed, the smoothed predicted density values ​​are integrated into a complete time-space three-dimensional tensor to form a predicted density field in an undisturbed state as the final output of step S3. The processing method is to combine the smoothed predicted density values ​​of each time point and each spatial grid unit into a three-dimensional data structure, covering all time points and spatial positions. Specifically, the smoothed predicted density values ​​of all time points and spatial grid units are uniformly stored and organized to ensure the integrity of the data structure and generate a predicted density field.

[0084] The predicted density field, as the benchmark data of step S4, needs to be provided in the form of a three-dimensional tensor so that it is consistent with the observed density field generated in step S1 and is convenient for point-by-point comparison. The unified three-dimensional data structure can directly support subsequent analysis, ensure the consistency of the data format and the smoothness of processing, and provide efficient and reliable input for difference analysis.

[0085] Step S3 reconstructs the natural evolution trend using low-frequency gradient components and performs smoothing to generate a predicted density field without interference. This predicted density field accurately reflects the natural evolution of the density field over time without human intervention, providing reliable benchmark data for step S4, which is used to perform difference analysis with the observed density field and locate the human interference area.

[0086] Step S1 uses multi-phase UAV remote sensing point cloud data to generate a spatiotemporal continuous density field, laying the foundation for subsequent analysis; Step S2 decomposes the change of the density field into a low-frequency gradient component of natural evolution and a high-frequency gradient component of artificial mutation through time axis gradient analysis; Step S3 generates a predicted density field under interference-free conditions based on the natural evolution component, providing a benchmark for natural changes. However, it is difficult to accurately identify the human interference area by simply comparing the predicted density field with the actual observed density field, especially in the face of complex surface changes and spatial heterogeneity. Therefore, step S4 introduces methods such as multi-scale difference calculation, comprehensive difference fusion, and dynamic threshold setting. By analyzing the difference between the predicted and observed density fields at different spatial scales, and combining the regional historical change characteristics with the spatial environment characteristics to dynamically adjust the threshold, the precise positioning of the human interference area can be achieved. The output of step S4 will provide reliable human interference area identification results for the subsequent step S5, paving the way for further causal analysis.

[0087] Step S4 includes the following contents:

[0088] S4.1, Multi-scale difference analysis:

[0089] In the process of locating the human interference area, it is first necessary to calculate the difference between the observed density field and the predicted density field at different spatial scales to capture the different levels of changes from local details to global trends. The processing method is to select a series of spatial scales, each of which corresponds to a Gaussian smoothing parameter, and then perform Gaussian smoothing on the observed density field and the predicted density field at each spatial scale to obtain the smoothed observed density field and the smoothed predicted density field. Then, for each time point and each spatial grid unit, the absolute value of the difference between the smoothed observed density value and the smoothed predicted density value at the current spatial scale is calculated, and recorded as the difference value at the time point and spatial grid unit at the spatial scale.

[0090] Differences at different spatial scales can reflect different range characteristics of changes. For example, small spatial scales are suitable for capturing local mutation characteristics, while large spatial scales are suitable for reflecting overall trend changes. Multi-scale analysis improves the comprehensiveness of difference detection and adaptability to complex changes by covering multiple spatial ranges. By integrating information at multiple spatial scales, the diverse manifestations of human disturbance areas can be more accurately identified.

[0091] S4.2, weighted fusion of inter-scale consistency:

[0092] After completing the multi-scale difference calculation, the difference fields at different spatial scales need to be fused into a comprehensive difference field to comprehensively evaluate the significance of the change. This is done by using consistency weighted fusion. The specific steps are as follows: First, for the difference field at each spatial scale, calculate its entropy value in the local neighborhood. The smaller the entropy value, the higher the consistency of the difference field at that spatial scale in the local neighborhood; then, calculate the weight of each spatial scale based on the entropy value. The smaller the entropy value, the larger the corresponding weight; finally, perform weighted average of the difference values ​​at each spatial scale according to the corresponding weights to obtain the comprehensive difference value of each time point and spatial grid unit.

[0093] Consistency weighted fusion can automatically enhance the contribution of spatial scales with consistent differences in local neighborhoods, while reducing the interference of spatial scales with noise or low consistency. The comprehensive difference field has higher robustness when fusing multi-scale information, and can more accurately reflect the true intensity of the change between the actual observed density field and the predicted density field, providing a reliable basis for the subsequent judgment of the human interference area.

[0094] S4.3, Dynamic Threshold Setting:

[0095] In order to adapt to the historical change characteristics and spatial heterogeneity of different regions, a dynamic threshold is set for each spatial grid unit as a basis for judging human interference. The dynamic threshold is calculated by combining the historical change amplitude and spatial heterogeneity. The specific steps are: first, the historical change amplitude of each spatial grid unit is calculated, that is, the average absolute difference of the observed density value change of the spatial grid unit at the historical time point; second, the spatial heterogeneity of the spatial grid unit is calculated, and the local Moran's I index is used to express it, reflecting the degree of spatial difference between the spatial grid unit and its neighborhood; finally, the historical change amplitude and spatial heterogeneity are weighted and summed by adjusting parameters to generate the dynamic threshold of each spatial grid unit.

[0096] The local Moran's I index is a statistical indicator used to measure spatial autocorrelation, which aims to analyze the similarity or difference between a specific spatial unit and its neighboring spatial units in a certain attribute value. Specifically, the local Moran's I index is used to quantify the spatial correlation between the observed density value of each spatial grid unit and the observed density values ​​of other spatial grid units in its neighborhood, thereby reflecting the spatial heterogeneity of the spatial grid unit.

[0097] Calculation principle and process:

[0098] The calculation of the local Moran's I index is based on the spatial weight matrix and the deviation analysis of the attribute values. For a specific spatial grid cell, the calculation process includes the following steps:

[0099] Determine the neighborhood range: Select the neighborhood of the target spatial grid cell, usually defined by a fixed distance or adjacent cells (such as a 5x5 grid), to form a set of spatial grid cells within the neighborhood.

[0100] Calculate attribute deviation: Calculate the difference between the observed density value of the target spatial grid cell and the average value of the observed density values ​​of all spatial grid cells in the neighborhood, and at the same time calculate the difference between the observed density value of each spatial grid cell in the neighborhood and the same average value.

[0101] Introduce spatial weights: assign weights based on the spatial relationship between the target spatial grid cell and each spatial grid cell in the neighborhood. Usually, the weight is inversely proportional to the distance or is set based on the adjacency relationship (such as adjacent is 1 and non-adjacent is 0).

[0102] Calculate the index: multiply the observed density deviation of the target spatial grid cell by the sum of the weighted deviations of each spatial grid cell in the neighborhood, and then divide it by the sum of the squares of the observed density deviations in the neighborhood to obtain the local Moran's I index.

[0103] The value of the local Moran's I index can be positive, negative, or close to zero:

[0104] Positive value: indicates that the observed density value of the target spatial grid cell shows a similar trend to the observed density values ​​of the spatial grid cells in its neighborhood, that is, the spatial aggregation is strong (such as high value aggregation or low value aggregation).

[0105] Negative values: indicate that the observed density value of the target spatial grid cell shows an opposite trend to the observed density values ​​of the spatial grid cells in its neighborhood, that is, the spatial heterogeneity is high (e.g., high values ​​are adjacent to low values).

[0106] Close to zero: Indicates that there is no significant spatial correlation between the observed density values ​​of the target spatial grid cell and the spatial grid cells in the neighborhood.

[0107] The local Moran's I index is used to measure the spatial heterogeneity of each spatial grid cell. Specifically, by calculating this index, it can be determined whether the observed density value of a spatial grid cell is significantly different from the observed density value of the spatial grid cells in its neighborhood. Areas with high spatial heterogeneity (such as the junction of vegetation and bare land) usually require a higher judgment threshold to avoid misjudging natural spatial changes as human interference; while areas with low spatial heterogeneity (such as uniform vegetation coverage areas) can use lower thresholds to improve detection sensitivity. This method ensures the adaptability of dynamic thresholds and enhances the accuracy of locating human interference areas.

[0108] The natural change characteristics and spatial environment of different regions are different. Fixed thresholds cannot adapt to this diversity. Dynamic thresholds are adaptively adjusted according to the specific characteristics of each spatial grid unit, which can improve the accuracy of judgment. Dynamic threshold setting reduces the possibility of misjudgment, especially in areas with frequent changes or high spatial heterogeneity, to avoid misidentifying natural changes as human interference.

[0109] S4.4, Locate the area of ​​human interference:

[0110] After generating the comprehensive difference field and dynamic threshold, the human interference area is located by comparing the comprehensive difference field with the dynamic threshold. For each time point and each spatial grid unit, it is determined whether its comprehensive difference value is greater than the dynamic threshold of the corresponding spatial grid unit. If it is greater than the dynamic threshold, it is determined that there is human interference at this time point and spatial grid unit, and the interference flag is recorded as 1; if it is less than or equal to the dynamic threshold, the interference flag is recorded as 0. Finally, the interference flags of all time points and spatial grid units are integrated to generate a human interference area map.

[0111] The comprehensive difference value reflects the degree of deviation between the observed density field and the predicted density field, and the dynamic threshold provides a judgment standard based on regional characteristics. The combination of the two can accurately identify abnormal changes caused by human intervention. Through quantitative comparison, the automatic positioning of human interference areas is achieved, which improves the practicality and processing efficiency of the monitoring method.

[0112] Step S4 generates a map of human interference areas through the process of multi-scale difference calculation, weighted fusion of inter-scale consistency, dynamic threshold setting and positioning of human interference areas. The map of human interference areas accurately locates the human interference areas based on the significant difference between the observed density field and the predicted density field, providing key input for subsequent causal analysis. This step uses multi-scale analysis and dynamic threshold setting methods to enhance the robustness of human interference detection and its adaptability to complex environments, ensuring the scientificity and practicality of the monitoring method.

[0113] Step S4 locates the human disturbance area and generates a human disturbance area map by multi-scale difference analysis of the observed density field and the predicted density field, combined with dynamic thresholds. The task of step S5 is to use the disturbance area map output by step S4, combined with geospatial data and historical activity information, to construct a causal reasoning network, analyze the causal sources and propagation paths of the human disturbance area, and provide decision support for natural resource management.

[0114] Step S5 includes the following contents:

[0115] S5.1, Data preparation and integration:

[0116] In the data preparation and integration stage, three types of input data are first collected and processed: maps of human disturbance areas, geospatial data, and historical activity information. The maps of human disturbance areas are generated by step S4 and are used to indicate whether there is human disturbance at a specific time and space grid coordinate. The disturbance flag is 1 for the presence of disturbance and 0 for the absence of disturbance. The geospatial data include the topographic features, land use types, vegetation coverage, and water system distribution information of the study area, providing a basis for the spatial association of the subsequent causal network. The historical activity information records the time, location, and intensity of human activities in the study area, such as mining, logging, agricultural development, and infrastructure construction, and is used to trace the source of disturbance. The processing process aligns these three types of data within the time series and spatial grid framework to ensure that the time points and spatial coordinates of all data are consistent for subsequent analysis.

[0117] Data integration provides a unified basis for causal analysis, ensuring accurate correspondence in time and space dimensions and avoiding analytical biases caused by inconsistent data. Comprehensive incorporation of geospatial data and historical activity information can more realistically reflect the external driving factors of interference and improve the scientificity and practicality of analysis.

[0118] S5.2, Causal network construction:

[0119] After completing data preparation and integration, a causal network is constructed to reflect the interference sources, interference areas and the causal relationship between them. The nodes in the causal network are divided into interference source nodes and interference area nodes. The interference source node is defined based on historical activity information. Each human activity record corresponds to an interference source node, which contains the activity type, occurrence time and spatial location. The interference area node is defined based on the human interference area map. The grid unit marked as interference at each time point is defined as the interference area node. The edges in the causal network include temporal causal edges, spatial causal edges and source-area edges. The temporal causal edge connects the interference area nodes of the same grid at adjacent time points, indicating the propagation of interference in time; the spatial causal edge connects the spatially adjacent interference area nodes at the same time point, indicating the diffusion of interference in space; the source-area edge connects the interference source node and the interference area node, indicating the direct impact of human activities on a specific area. Finally, a directed graph containing all nodes and edges is constructed to provide a structural basis for causal reasoning. The reason for the technical features is that by clearly defining nodes and edges, the causal network can intuitively represent the spatiotemporal propagation and causal sources of interference. The causal network structure provides a visual and systematic framework for causal reasoning, which helps to identify complex causal chains.

[0120] Causal Reasoning:

[0121] After the causal network is constructed, causal reasoning is performed to evaluate and screen the edges in the causal network. Causal reasoning is divided into three parts: temporal causal reasoning, spatial causal reasoning, and source-region causal reasoning.

[0122] Temporal causal reasoning: The Pearson correlation coefficient is used to evaluate the temporal correlation between nodes in the interference area at adjacent time points. The calculation method is to first calculate the covariance of the interference signs at two time points, then calculate the standard deviation of the interference signs at two time points respectively, and then divide the covariance by the product of the two standard deviations to obtain the correlation coefficient to measure the continuity of the interference in time. If the calculated correlation coefficient is greater than the preset time threshold, the corresponding temporal causal edge is retained.

[0123] Spatial causal reasoning: The Moran's I index is used to evaluate the spatial correlation between adjacent grids at the same time point. The calculation method is to first calculate the deviation of the interference mark between the target grid and the neighboring grid, multiply these deviations two by two and sum them, and then standardize them with the total deviation of the interference mark of all grids in the study area to obtain the Moran's I index, which measures the degree of spatial aggregation or diffusion of interference. If the calculated Moran's I index is greater than the preset spatial threshold, the corresponding spatial causal edge is retained.

[0124] Source-region causal reasoning: Granger causality test is used to evaluate the causal impact of interference source nodes on interference region nodes. The calculation method is to first build a time series model that includes interference source activities and a time series model that does not include interference source activities, compare the goodness of fit of the two models, and calculate the p-value under the significance level. If the p-value is less than the preset significance level, the corresponding source-region edge is retained.

[0125] The use of mature statistical methods such as Pearson correlation coefficient, Moran's I index and Granger causality test for causal reasoning ensures the scientificity and objectivity of the analysis. These methods can quantitatively evaluate the strength and significance of causal relationships and provide a reliable basis for screening edges in the causal network.

[0126] S5.3, Causal Path Analysis:

[0127] After completing the causal reasoning, causal path analysis is performed to identify the impact path of the interference source on the interference area. The processing method is to calculate the shortest path from all interference source nodes to the interference area node for each interference area node in the constructed and screened causal network, and the path weight is defined as the inverse of the association strength. The calculation process is to start from the interference source node, traverse along the directed edges step by step, accumulate the weights on the path, and find the path with the smallest total weight as the shortest path. Based on the shortest path, identify the interference source and its propagation path that has the greatest impact on the interference area node. Finally, the causal sources of all interference area nodes are counted, and the main interference sources and their impact ranges are aggregated to identify.

[0128] The shortest path algorithm can efficiently identify the optimal causal chain between the interference source and the interference area, revealing the propagation mechanism of the interference. Through causal path analysis, the source and propagation process of the interference can be intuitively displayed, providing targeted intervention suggestions for resource management.

[0129] Step S5 constructs a causal reasoning network by integrating the human disturbance area map, geospatial data and historical activity information generated in step S4, and realizes the analysis of the causal source and propagation path of the human disturbance area. The network clearly reveals the impact of human activities as a source of interference on natural resources and its spatiotemporal propagation law. The main interference sources and their impact range are identified through causal path analysis, providing a targeted basis for protection and intervention measures in dynamic monitoring of natural resources.

[0130] The above formulas are all dimensionless and numerical calculations. The formula is a formula for the most recent real situation obtained by collecting a large amount of data and performing software simulation. The preset parameters in the formula are set by technicians in this field according to actual conditions.

[0131] It should be noted that the system of the present invention can be deployed on the device itself to realize embedded applications, and can also be run on a PC or other terminal with a user interface, thereby meeting a variety of hardware environments and usage requirements.

[0132] The above description is only by way of illustration of certain exemplary embodiments of the present invention. It is undoubted that those skilled in the art can modify the described embodiments in various ways without departing from the spirit and scope of the present invention. Therefore, the above drawings and descriptions are illustrative in nature and should not be construed as limiting the scope of protection of the claims of the present invention.

[0133] It should be noted that, in this article, if there are relational terms such as first and second, etc., they are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, article or device. In the absence of further restrictions, the elements defined by the sentence "including a..." do not exclude the existence of other identical elements in the process, method, article or device including the elements.

[0134] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art who is familiar with the present technical field can easily think of changes or substitutions within the technical scope disclosed in the present application, which should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be based on the protection scope of the claims.

Claims

1. A UAV remote sensing Gaussian splash change detection method for dynamic monitoring of natural resources, characterized in that: Includes steps: S1. Obtain multi-temporal UAV remote sensing point cloud data and convert it into a spatiotemporal continuous observation density field; S2. Perform time axis gradient analysis on the observed density field, separate the low frequency gradient components of natural evolution and the high frequency gradient components of artificial mutation, and generate a classification map of change patterns; S3. Based on the naturally evolving low-frequency gradient component, obtain the predicted density field under the undisturbed state; S4. Calculate the difference between the observed density field and the predicted density field at different spatial scales, and fuse the differences between scales into comprehensive differences through consistency weighted fusion. Combine the historical change amplitude of the region and the spatial heterogeneity to dynamically set the judgment threshold, and locate the human interference area based on the comparison between the comprehensive difference and the dynamic threshold. S5. Combine geospatial data and historical activity information to construct a causal reasoning network to analyze the causal sources and propagation paths of human disturbance areas.

2. The method for detecting changes in Gaussian splashes using remote sensing by unmanned aerial vehicles for dynamic monitoring of natural resources according to claim 1 is characterized in that: Step S1 includes the following contents: The multi-temporal UAV remote sensing point cloud data are denoised, and a density-based clustering method is used to remove outliers to generate denoised point cloud data. All denoised point cloud data are aligned to a unified geographic coordinate system through a registration method. The study area is divided into a uniform two-dimensional grid and subdivided into multiple altitude layers in the vertical direction. The denoised point cloud data are assigned to the corresponding two-dimensional grid cells and altitude layers according to the three-dimensional coordinates. The local observation density value is calculated for the point set in each altitude layer using the Gaussian kernel density estimation method, and the comprehensive observation density value is generated by weighted summation. The comprehensive observation density value is linearly interpolated in the time dimension to fill the missing phase data, and Gaussian filtering is applied for spatiotemporal smoothing to finally generate the observation density field.

3. The method for detecting changes in Gaussian splashes using remote sensing by unmanned aerial vehicles for dynamic monitoring of natural resources according to claim 2 is characterized in that: Step S2 includes the following contents: The observed density field data is subjected to time axis gradient analysis. The density change gradient is calculated by the central difference method for the observed density value sequence of each spatial grid unit. The density change gradient sequence is converted into frequency domain spectrum data by discrete Fourier transform. The spectrum data is separated into low-frequency components and high-frequency components according to a preset frequency threshold. The low-frequency components and high-frequency components are respectively reconstructed into low-frequency gradient components and high-frequency gradient components in the time domain by inverse discrete Fourier transform. The relative intensity ratio of the high-frequency gradient component to the low-frequency gradient component is calculated and compared with the classification threshold to generate a change pattern classification map. Finally, the low-frequency gradient component field, high-frequency gradient component field and change pattern classification map are output as input for predicting the density field and locating the human interference area.

4. The method for detecting changes in Gaussian splashes using remote sensing by unmanned aerial vehicles for dynamic monitoring of natural resources according to claim 3 is characterized in that: Step S3 includes the following contents: Based on the separated low-frequency gradient components, the observed density values ​​at the initial time points provided by the observed density field are used to reconstruct the predicted density values ​​at subsequent time points by recursively accumulating the product of the low-frequency gradient components and the time interval. The reconstructed predicted density field is smoothed by applying Gaussian filtering in the time dimension. The smoothed predicted density values ​​are generated by taking Gaussian weighted average of the predicted density values ​​at adjacent time points in the time window. The smoothed predicted density values ​​are integrated into a time-space three-dimensional tensor to form a predicted density field in an undisturbed state.

5. The method for detecting Gaussian splash changes using remote sensing by unmanned aerial vehicles for dynamic monitoring of natural resources according to claim 4 is characterized in that: Step S4 includes the following contents: For the observed density field and the predicted density field, the absolute value of the difference between the smoothed observed density value and the smoothed predicted density value is calculated by Gaussian smoothing at multiple spatial scales as the difference value. The consistent weighted fusion method is used to calculate the entropy value of the difference field in the local neighborhood at each spatial scale, and the weight is determined according to the entropy value. The comprehensive difference value is generated by weighted average. The dynamic threshold is set based on the spatial heterogeneity calculated by combining the average absolute difference of the historical observed density value changes of each spatial grid unit and the local Moran's I index. The map of the human interference area is generated by comparing the comprehensive difference value with the dynamic threshold.

6. The method for detecting changes in Gaussian splashes using remote sensing by unmanned aerial vehicles for dynamic monitoring of natural resources according to claim 5 is characterized in that: Step S5 includes the following contents: The human interference area map, geospatial data and historical activity information are integrated and aligned in the time and space dimensions. After defining the interference source nodes based on the historical activity information and the interference area nodes based on the human interference area map, a causal network containing temporal causal edges, spatial causal edges and source-area edges is constructed. Based on the screened causal network, the shortest path algorithm is used to identify the interference sources and their propagation paths, and the main interference sources and the impact range are counted.

7. The method for detecting changes in Gaussian splashes using remote sensing by unmanned aerial vehicles for dynamic monitoring of natural resources according to claim 6 is characterized in that: Step S5 also includes the following contents: The temporal correlation of nodes in the interference region at adjacent time points is calculated using the Pearson correlation coefficient and compared with the time threshold to screen the temporal causal edges.

8. The method for detecting changes in Gaussian splashes using remote sensing by unmanned aerial vehicles for dynamic monitoring of natural resources according to claim 6, characterized in that: Step S5 also Includes the following: The spatial correlation between adjacent grids at the same time point is calculated by the Moran's I index and compared with the spatial threshold to screen the spatial causal edges.

9. The method for detecting changes in Gaussian splashes using remote sensing by unmanned aerial vehicles for dynamic monitoring of natural resources according to claim 6, characterized in that: Step S5 also Includes the following: The causal impact of the interference source node on the interference region node is calculated through the Granger causality test and compared with the significance level to screen the source-region edge.

Images