UAV Remote Sensing Gaussian Splashing Change Detection Method for Dynamic Monitoring of Natural Resources

Abstract

The present invention discloses a UAV remote sensing Gaussian splash change detection method for dynamic monitoring of natural resources, which specifically relates to the field of natural resource monitoring and is used to solve the problem that it is difficult to distinguish natural evolution from human interference. It generates a spatiotemporally continuous observation density field through the conversion of multi-temporal point cloud data, combines time-axis gradient analysis to separate natural and human change features, uses low-frequency gradient components to predict the density field in the interference-free state, and accurately locates human interference areas through multi-scale difference fusion and dynamic threshold setting. Finally, a causal inference network is constructed to reveal the source and propagation path of interference. At the same time, by integrating geospatial data with historical activity information, the causal association analysis of interference events is realized, providing targeted basis for ecological protection decision-making, effectively improving the timeliness of natural resource supervision, and providing efficient support for resource management and environmental protection.

Description

Technical Field

[0001] The present invention relates to the field of natural resource monitoring. More specifically, the present invention relates to an unmanned aerial vehicle remote sensing Gaussian splash change detection method for natural resource dynamic monitoring. Background Art

[0002] In the field of natural resource dynamic monitoring, the existing technology generally adopts a density difference detection method based on a fixed threshold. Its core limitation lies in the inability to effectively model the essential differences between natural evolution processes and human intervention behaviors. Natural evolution usually shows a progressive change that conforms to ecological laws on a long time scale, while human intervention often triggers local mutations that violate natural laws in a short time. Although they may show approximate characteristics in the change amplitude of the density field, their dynamic patterns in the time dimension are completely different. Due to the lack of the ability to quantitatively model the rate of temporal change, the existing methods misclassify slow natural processes and sudden human events as the same type of change; 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 scenarios such as seasonal withering and flourishing of vegetation and natural river course changes, resulting in the dual problems of increased false alarm rates and missed detection of key events. More seriously, the existing technology cannot analyze the causal relationships between change events. For example, it cannot identify whether the vegetation degradation in a certain area is indirectly caused by adjacent mining activities, resulting in subsequent treatment measures lacking pertinence. To solve the above problems, a technical solution is provided now. Summary of the Invention

[0003] To overcome the above-mentioned defects of the existing technology, an embodiment of the present invention provides an unmanned aerial vehicle remote sensing Gaussian splash change detection method for natural resource dynamic monitoring. By converting multi-temporal point cloud data, a spatio-temporally continuous observed density field is generated. Combining time-axis gradient analysis to separate natural and human-induced change features, using low-frequency gradient components to predict the density field in a non-interference state, and accurately locating human interference areas through multi-scale difference fusion and dynamic threshold setting. Finally, a causal inference network is constructed to reveal the source and propagation path of interference; at the same time, by integrating geospatial data and historical activity information, the causal relationship analysis of interference events is realized, providing a targeted basis for ecological protection decision-making, providing efficient support for resource management and environmental protection, and solving the problems raised in the above background art.

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

[0005] S1. Obtain multi-temporal unmanned aerial vehicle remote sensing point cloud data and convert it into a spatio-temporally continuous observed density field;

[0006] S2. Conduct 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 human-induced mutations, and generate a change pattern classification map;

[0007] S3. Obtain the predicted density field in an interference-free state based on the low-frequency gradient component of natural evolution;

[0008] S4. Calculate the difference between the observed density field and the predicted density field at different spatial scales, fuse them into a comprehensive difference through the consistency weighting of the differences between scales, dynamically set the judgment threshold in combination with the regional historical change amplitude and spatial heterogeneity, and locate the human interference area according to the comparison between the comprehensive difference and the dynamic threshold;

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

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

[0011] Perform denoising processing on multi-temporal UAV remote sensing point cloud data, use density-based clustering method to remove outliers to generate denoised point cloud data, align all denoised point cloud data to a unified geographic coordinate system through registration method, divide the study area into uniform two-dimensional grids and subdivide them into multiple height layers in the vertical direction, assign the denoised point cloud data to the corresponding two-dimensional grid cells and height layers according to the three-dimensional coordinates, calculate the local observed density value for the point set in each height layer using the Gaussian kernel density estimation method and generate the comprehensive observed density value through weighted summation, perform linear interpolation on the comprehensive observed density value in the time dimension to fill in the missing temporal data and apply Gaussian filtering for spatio-temporal smoothing processing, and finally generate the observed density field.

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

[0013] Perform time-axis gradient analysis on the observed density field data, calculate the density change gradient for the sequence of observed density values of each spatial grid cell using the central difference method, perform discrete Fourier transform on the density change gradient sequence to convert it into frequency-domain spectrum data, separate the spectrum data into low-frequency components and high-frequency components according to the preset frequency threshold, perform inverse discrete Fourier transform on the low-frequency components and high-frequency components respectively to reconstruct them into low-frequency gradient components and high-frequency gradient components in the time domain, calculate the relative intensity ratio of the high-frequency gradient components to the low-frequency gradient components and compare it with the classification threshold to generate a change pattern classification map, and finally output the low-frequency gradient component field, high-frequency gradient component field and change pattern classification map as the input for predicting the density field and locating human interference areas.

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

[0015] Based on the separated low-frequency gradient components, using the observed density values at the initial time point provided by the observed density field, the predicted density values at subsequent time points are reconstructed 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 performing Gaussian weighted averaging on the predicted density values at adjacent time points within a time window. The smoothed predicted density values are integrated into a three-dimensional time-space tensor to form a predicted density field in a non-interference state.

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

[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 as the difference value by Gaussian smoothing at multiple spatial scales. The entropy value of the difference field within the local neighborhood at each spatial scale is calculated using the consistency weighted fusion method, and the weights are determined based on the entropy value and then weighted averaged to generate a comprehensive difference value. A dynamic threshold is set by combining the average absolute difference in the historical observed density value changes of each spatial grid cell and the spatial heterogeneity calculated by the local Moran's I index. The artificial interference area map is generated by comparing the comprehensive difference value with the dynamic threshold.

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

[0019] Integrate the artificial interference area map, geospatial data, and historical activity information and align them in the spatio-temporal dimension. After defining the interference source nodes based on the historical activity information and the interference area nodes based on the artificial interference area map, construct a causal network including time causal edges, spatial causal edges, and source-region edges. Based on the filtered causal network, use the shortest path algorithm to identify the interference sources and their propagation paths and count the main interference sources and the affected ranges.

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

[0021] Calculate the temporal correlation degree of the interference area nodes at adjacent time points through the Pearson correlation coefficient and compare it with the time threshold to screen the time causal edges.

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

[0023] Calculate the spatial correlation degree between adjacent grids at the same time point through the Moran's I index and compare it with the spatial threshold to screen the spatial causal edges.

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

[0025] Calculate the causal influence of interference source nodes on nodes in the interference area through Granger causality test and compare it with the significance level to screen source-region edges.

[0026] Technical effects and advantages of the UAV remote sensing Gaussian splash change detection method for natural resource dynamic monitoring of the present invention:

[0027] Through the UAV remote sensing Gaussian splash change detection method for natural resource dynamic monitoring, the present invention significantly improves the accuracy and scientific nature of natural resource change monitoring. Aiming at the limitations of the fixed threshold method in the prior art, such as being difficult to distinguish natural evolution from human intervention, being vulnerable to noise interference, and lacking causal analysis, the present invention generates a spatio-temporally continuous observation density field through multi-temporal point cloud data conversion, combines time-axis gradient analysis to separate natural and human change features, uses low-frequency gradient components to predict the density field in the interference-free state, and accurately locates human interference areas through multi-scale difference fusion and dynamic threshold setting. Finally, a causal inference network is constructed to reveal the source and propagation path of interference. This technical solution overcomes the deficiencies of traditional methods in modeling the temporal change rate and confusing pseudo-change signals, reduces the false alarm rate and missed detection rate, and enhances the ability to distinguish complex scenes. At the same time, by integrating geospatial data and historical activity information, the present invention realizes the causal association analysis of interference events and provides a targeted basis for ecological protection decision-making. Brief Description of the Drawings

[0028] Figure 1 It is a flow chart of the UAV remote sensing Gaussian splash change detection method for natural resource dynamic monitoring of the present invention. Detailed Embodiments

[0029] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0030] Embodiment 1: Figure 1 The UAV remote sensing Gaussian splash change detection method for natural resource dynamic monitoring of the present invention is given, including:

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

[0032] S2. Conduct time-axis gradient analysis on the observation density field, separate the low-frequency gradient components of natural evolution and the high-frequency gradient components of human mutations, and generate a change pattern classification map.

[0033] S3. Obtain the predicted density field in the interference-free state based on the low-frequency gradient component of natural evolution.

[0034] S4. Calculate the difference between the observed density field and the predicted density field at different spatial scales, and through the consistency weighted fusion of the differences between scales into a comprehensive difference, dynamically set the judgment threshold in combination with the regional historical change amplitude and spatial heterogeneity, and locate the human interference area according to the comparison between the comprehensive difference and the dynamic threshold.

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

[0036] In the field of dynamic monitoring of natural resources, it is crucial to timely grasp phenomena such as vegetation cover change, river morphology adjustment, and land use pattern transformation for ecological protection and resource management. These changes are often jointly affected by the natural evolution process and human intervention behaviors. Therefore, it is necessary to obtain high-resolution multi-temporal point cloud data through unmanned aerial vehicle (UAV) remote sensing technology to construct an observed density field data that can reflect the spatio-temporal dynamics of surface features, so as to provide a data basis for subsequent distinguishing natural and human changes, locating interference areas, and analyzing causal relationships. And step S1 is precisely through the processing of multi-temporal UAV remote sensing point cloud data to generate spatio-temporally continuous observed density field data, laying a 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 a set of point clouds collected at different time points. Each group of point cloud data contains the three-dimensional coordinate information of a large number of points, that is, the abscissa in the horizontal direction, the ordinate in the horizontal direction, and the elevation coordinate in the vertical direction. These multi-temporal UAV remote sensing point cloud data are the basic input data for subsequent conversion of the observed density field.

[0040] Multi-temporal UAV remote sensing point cloud data refers to a set of point cloud data collected at multiple different time points by remote sensing equipment carried by UAVs. Specifically, the remote sensing equipment (such as lidar or high-resolution cameras) repeats the acquisition of the three-dimensional spatial information of the surface of the same study area at different time periods (such as daily, monthly, or annually) during the flight of the UAV. These information are recorded in the form of point clouds, and each point contains the abscissa in the horizontal direction, the ordinate in the horizontal direction, and the elevation coordinate in the vertical direction, and may also include additional attributes, such as reflection intensity or color values. The so-called multi-temporal emphasizes that the acquisition of data has the characteristics of a time series and can reflect the dynamic process of surface features changing over time, such as phenomena such as vegetation growth, river course change, or human development.

[0041] S1.2, Point cloud preprocessing:

[0042] Denoise each set of multi-temporal UAV remote sensing point cloud data. Specifically, use the density-based clustering method. By analyzing the distance between points, identify and remove outliers that deviate from the main point cloud distribution, and generate denoised point cloud data. After denoising, use the registration method. By comparing the spatial position differences between the point cloud data collected at different time points, align all denoised point cloud data into a unified geographic coordinate system to ensure the consistency of the multi-temporal point cloud data in spatial position.

[0043] S1.3, Spatial grid division:

[0044] Divide the study area into uniform two-dimensional grids. The horizontal and vertical spacings of each grid cell are fixed values set in advance, used to represent specific spatial positions. On this basis, along the vertical direction, i.e., the elevation direction, each grid cell is further divided into multiple height layers, and the interval between adjacent height layers is a fixed value set in advance, forming a three-dimensional spatial division structure.

[0045] The two-dimensional grid division discretizes the continuous space of the study area into computable units, facilitating subsequent density calculation and spatial distribution analysis. The height layer subdivision can reflect the distribution differences of the point cloud data in the vertical direction, such as distinguishing the vegetation canopy from surface features, enhancing the comprehensiveness of spatial description.

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

[0047] For each set of denoised multi-temporal UAV remote sensing point cloud data, according to the three-dimensional coordinates of the points, allocate all points to the corresponding two-dimensional grid cells and their height layers. For the point set within a specific height layer of each grid cell, use the Gaussian kernel density estimation method to calculate the local density. The specific calculation process is as follows: First, determine the distance between each point in this height layer and the geometric center of the grid cell, and then assign a weight to each point according to the distance through the Gaussian function, where the weight size is controlled by a preset smoothing bandwidth parameter; then, sum up the weight values of all points in this height layer and perform normalization processing to obtain the local observation density value of this height layer. After that, perform weighted summation on the local observation density values of each grid cell in all height layers to generate the comprehensive observation density value of this grid cell at the current time point, where the weighting coefficient is adjusted in advance according to the land cover type (such as vegetation or buildings).

[0048] The Gaussian kernel density estimation method reduces local noise in the point cloud data distribution through smoothing processing, providing a continuous and stable density representation. The weighted integration of the local density of the height layers comprehensively considers the distribution characteristics of the land cover in the vertical direction, enhancing the representation ability of the comprehensive observation density value for surface features.

[0049] S1.5, Spatiotemporal Continuity Processing:

[0050] For the missing time points in the multi-temporal UAV remote sensing point cloud data, the comprehensive observation density values of the adjacent front and rear time points are used, and the comprehensive observation density values of the missing time points are calculated by the linear interpolation method to fill the time series gap. After interpolation, the Gaussian filtering method is applied to the comprehensive observation density values of all time points and spatial positions for smoothing processing in both the time dimension and the spatial dimension to generate the smoothed observation density field data.

[0051] The linear interpolation method calculates the missing values by inferring the trend of adjacent data, ensuring the integrity of the time series and supporting subsequent time dimension-based analysis. The smoothing processing of the Gaussian filter in the spatio-temporal dimension reduces the influence of short-term fluctuations and spatial noise, improving the stability and continuity of the observation density field data.

[0052] S1.6, Observation Density Field Output:

[0053] Finally, a spatio-temporally continuous observation density field data structure is generated. This observation density field contains the smoothed comprehensive observation density values of all time points and all grid cells, forming a three-dimensional data set with dimensions of time, number of grid rows, and number of grid columns for subsequent analysis.

[0054] Step S1 has generated spatio-temporally continuous observation density field data through the processing of multi-temporal UAV remote sensing point cloud data, comprehensively reflecting the distribution characteristics of surface features in the time and space dimensions. Step S2 takes over this observation density field data, aiming to separate the low-frequency gradient component of natural evolution and the high-frequency gradient component of human-induced mutations through time-axis gradient analysis and generate a change pattern classification map to reveal the internal laws and abnormal characteristics of natural resource changes.

[0055] Step S2 includes the following:

[0056] S2.1, Time-Axis Gradient Calculation:

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

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

[0059] S2.2, Frequency domain conversion:

[0060] Perform frequency domain conversion on the density change gradient sequence of each spatial grid cell. Using the discrete Fourier transform method, convert the density change gradient data in the time dimension into spectral data in the frequency dimension. This conversion process decomposes the density change gradient sequence into multiple sine wave components with different frequencies. Each sine 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 with different frequencies, facilitating the distinction of the periodic characteristics of the changes. Natural evolution usually shows slow changes at low frequencies, while human-induced sudden changes usually show rapid changes at high frequencies. Through this conversion, a clear frequency domain basis can be provided for the separation of low-frequency 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, according to the pre-set frequency threshold, divide the spectral data into low-frequency and high-frequency components. The determination of the frequency threshold is based on the time scale of the natural evolution of natural resources. By analyzing historical data, statistically obtain the maximum frequency value of natural changes. The spectral components below this frequency threshold are classified as low-frequency components, representing the characteristics of natural evolution; the spectral components above this frequency threshold are classified as high-frequency components, representing the characteristics of human-induced sudden changes.

[0064] By using the frequency threshold to achieve the separation of spectral components, it 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, enabling the characteristics of natural evolution and human-induced sudden changes to be clearly distinguished in the frequency domain.

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

[0066] Perform inverse discrete Fourier transform on the separated low-frequency and high-frequency components respectively, converting them from the frequency domain back to the time domain, generating the time series of the low-frequency gradient components and the time series of the high-frequency gradient components. This reconstruction process restores the changing characteristics of the low-frequency and high-frequency components in the time dimension, facilitating 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, enabling these components to 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 convenient for processing 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 cell at each time point, calculate the relative intensity ratio of the high-frequency gradient component to the low-frequency gradient component. 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 a pre-set classification threshold. If this relative intensity ratio is greater than the classification threshold, it is determined that the change pattern at this time point in this spatial grid cell is an artificial mutation; if this relative intensity ratio is less than or equal to the classification threshold, it is determined as natural evolution.

[0070] The relative intensity ratio can quantify the dominant degree 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, a low-frequency gradient component field, a high-frequency gradient component field, and a change pattern classification map are generated. The low-frequency gradient component field represents the natural evolution trend of each spatial grid cell in the time dimension. The high-frequency gradient component field and the change pattern classification map are used to identify and locate the areas of human intervention, reflecting the distribution and intensity of artificial mutations.

[0073] Step S2 receives the observed density field data generated by step S1. 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, it separates the low-frequency gradient components representing natural evolution and the high-frequency gradient components representing artificial mutations, and generates a change pattern classification map. These processing results provide key input data for step S3 to predict the observed density field in the interference-free state and step S4 to locate the areas of human intervention, ensuring the logic and integrity of the entire method.

[0074] Step S2 successfully separates the low-frequency gradient components of natural evolution and the high-frequency gradient components of artificial mutations 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 component of natural evolution and the high-frequency gradient component of artificial mutation by performing time-axis gradient analysis on the spatio-temporally continuous observed density field data, and generated a classification map of change patterns. The low-frequency gradient component reflects the natural change trend of the observed density field over time and is not affected by artificial interference. Taking this low-frequency gradient component as the input, step S3 generates a predicted density field in a non-interference state through reconstruction and smoothing, aiming to provide reference data for step S4 to perform difference analysis with the observed density field and locate the artificial interference areas.

[0076] Step S3 includes the following:

[0077] S3.1, Reconstruct the natural evolution trend using the low-frequency gradient component:

[0078] In the process of generating the predicted density field in a non-interference state, it is first necessary to reconstruct the natural evolution trend of the observed density field in the time dimension using the low-frequency gradient component. The low-frequency gradient component is obtained through time-axis gradient analysis in step S2, representing the rate of natural change of the observed density field over time, and the influence of artificial mutation has been removed. The processing method is to start from the observed density value at the initial time point and gradually accumulate the low-frequency gradient component to predict the observed density values at subsequent time points. The specific operation steps are as follows: 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 cell, calculate the predicted density value at the subsequent time point. The calculation method of the predicted density value is to add the product of the low-frequency gradient component and the time interval to the predicted density value at the previous time point, and recursively obtain 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 sudden change components of artificial interference. Therefore, the observed density field reconstructed based on the low-frequency gradient component can truly simulate the natural changes in a non-interference state. The recursive accumulation calculation method is simple and efficient, can continuously predict density values in the time series, and ensure the consistency of the prediction results with the natural evolution trend.

[0080] S3.2, Smoothing processing to improve stability:

[0081] After reconstructing the natural evolution trend using the low-frequency gradient component, since there may be minute noise or local discontinuities in the low-frequency gradient component, the directly obtained predicted density field needs to be further smoothed to enhance its spatio-temporal continuity and reliability. The processing method is to apply Gaussian filtering to the reconstructed predicted density field in the time dimension, generating a smoothed predicted density field by performing weighted averaging on the predicted density values at adjacent time points. Specifically: Select a fixed time window, and within this window, perform a weighted summation of the predicted density values at each time point. The weighting method uses the weights of a Gaussian distribution, and the weight values are determined according to the distance between the time point and the central point. The closer the distance, the higher the weight, and the sum of all weights is 1, finally obtaining the smoothed predicted density value.

[0082] S3.3, Generate the predicted density field:

[0083] After completing the smoothing process, integrate the smoothed predicted density values into a complete three-dimensional time-space tensor to form the predicted density field in a non-interference state, which is the final output of step S3. The processing method is to combine the smoothed predicted density values at each time point and each spatial grid cell into a three-dimensional data structure that covers all time points and spatial positions. Specifically: Uniformly store and organize the smoothed predicted density values for all time points and spatial grid cells to ensure the integrity of the data structure and generate the predicted density field.

[0084] The predicted density field, as the reference data for step S4, needs to be provided in the form of a three-dimensional tensor to be consistent with the observed density field generated in step S1, facilitating point-by-point comparison. The unified three-dimensional data structure can directly support subsequent analysis, ensuring the coherence of the data format and the smoothness of processing, providing efficient and reliable input for the difference analysis.

[0085] In step S3, by using the low-frequency gradient component to reconstruct the natural evolution trend and performing smoothing processing, the predicted density field in a non-interference state is finally generated. This predicted density field accurately reflects the natural evolution state of the density field over time without human intervention, providing reliable reference data for step S4 to perform difference analysis with the observed density field and locate the human interference areas.

[0086] Step S1 generates a spatio-temporally continuous density field using multi-temporal UAV remote sensing point cloud data, laying the foundation for subsequent analysis; Step S2 decomposes the changes in the density field into low-frequency gradient components of natural evolution and high-frequency gradient components of human-induced mutations through time-axis gradient analysis; Step S3 then generates a predicted density field in an undisturbed state based on the natural evolution component, providing a benchmark for natural changes. However, simply relying on the simple comparison between the predicted density field and the actual observed density field is difficult to accurately identify human-interfered areas, 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 differences between the predicted and observed density fields at different spatial scales and dynamically adjusting the threshold in combination with regional historical change characteristics and spatial environmental features, the precise location of human-interfered areas can be achieved. The output of Step S4 will provide a reliable identification result of human-interfered areas for the subsequent Step S5, paving the way for further causal analysis.

[0087] Step S4 includes the following:

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

[0089] In the process of locating human-interfered areas, it is first necessary to calculate the differences between the observed density field and the predicted density field at different spatial scales to capture different levels of changes from local details to global trends. The processing method is to select a series of spatial scales, each corresponding to a Gaussian smoothing parameter, and then perform Gaussian smoothing on the observed density field and the predicted density field respectively 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 cell, calculate the absolute value of the difference between the smoothed observed density value and the smoothed predicted density value at the current spatial scale, which is denoted as the difference value of this time point and spatial grid cell at this spatial scale.

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

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

[0092] After completing the multi-scale difference calculation, it is necessary to fuse the difference fields at different spatial scales into a comprehensive difference field to comprehensively evaluate the significance of changes. Consistency weighted fusion is adopted, and the specific steps are as follows: First, for the difference field at each spatial scale, calculate its entropy value within the local neighborhood. The smaller the entropy value, the higher the consistency of the difference field within the local neighborhood at that spatial scale. Then, calculate the weight of each spatial scale according to the entropy value. The smaller the entropy value, the larger the corresponding weight. Finally, perform weighted averaging on the difference values at each spatial scale according to the corresponding weights to obtain the comprehensive difference value for each time point and spatial grid cell.

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

[0094] S4.3, Dynamic threshold setting:

[0095] To adapt to the historical change characteristics and spatial heterogeneity of different regions, a dynamic threshold is set for each spatial grid cell as the basis for judging human interference. The dynamic threshold is calculated by combining the historical change amplitude and spatial heterogeneity. The specific steps are as follows: First, calculate the historical change amplitude of each spatial grid cell, that is, the average absolute difference in the observed density value of the spatial grid cell at historical time points. Second, calculate the spatial heterogeneity of the spatial grid cell, which is represented by the local Moran's I index and reflects the spatial difference degree between the spatial grid cell and its neighborhood. Finally, perform weighted summation on the historical change amplitude and spatial heterogeneity through adjustment parameters to generate the dynamic threshold for each spatial grid cell.

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

[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 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 the property deviation: Calculate the difference between the observed density value of the target spatial grid cell and the average of the observed density values of all spatial grid cells within the neighborhood. At the same time, calculate the difference between the observed density value of each spatial grid cell within the neighborhood and the same average.

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

[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 within the neighborhood, and then divide by the sum of the squares of the observed density deviations within 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 within its neighborhood, that is, strong spatial aggregation (such as high - value aggregation or low - value aggregation).

[0105] Negative value: Indicates 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 within its neighborhood, that is, high spatial heterogeneity (such as high - value adjacent to low - value).

[0106] Close to zero: Indicates that there is no significant spatial correlation between the observed density value of the target spatial grid cell and the observed density values of the spatial grid cells within 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 there is a significant difference between the observed density value of a certain spatial grid cell and the observed density values of the spatial grid cells within 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 uniformly vegetated areas) can use a lower threshold to improve detection sensitivity. This method ensures the adaptability of the dynamic threshold and enhances the accuracy of locating human - disturbed areas.

[0108] There are differences in the natural change characteristics and spatial environments of different regions. Using a fixed threshold cannot adapt to this diversity, while a dynamic threshold can adaptively adjust according to the specific characteristics of each spatial grid cell, which can improve the accuracy of judgment. The setting of the dynamic threshold reduces the possibility of misjudgment, especially in regions with frequent changes or high spatial heterogeneity, avoiding misidentifying natural changes as human interference.

[0109] S4.4, Locate the human interference areas:

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

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

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

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

[0114] Step S5 includes the following:

[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: the artificial interference area map, geospatial data, and historical activity information. The artificial interference area map is generated by step S4 and is used to indicate the presence or absence of artificial interference at specific time and spatial grid coordinates. The interference flag uses 1 to indicate the presence of interference and 0 to indicate no interference. The geospatial data includes the topographic features, land use types, vegetation cover, 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, etc., and is used to trace the source of interference. 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, facilitating subsequent analysis.

[0117] Data integration provides a unified basis for causal analysis, ensuring accurate correspondence in the time and space dimensions and avoiding analysis biases caused by inconsistent data. Comprehensively incorporating geospatial data and historical activity information can more realistically reflect the external driving factors of interference, improving the scientificity and practicality of the 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 relationships between them. The nodes in the causal network are divided into interference source nodes and interference area nodes. The interference source nodes are defined based on historical activity information, and each human activity record corresponds to an interference source node, including the activity type, occurrence time, and spatial location. The interference area nodes are defined based on the artificial interference area map, and each grid cell marked as interfered at each time point is defined as an interference area node. The edges in the causal network include time causal edges, spatial causal edges, and source-area edges. The time causal edges connect the interference area nodes of the same grid at adjacent time points, indicating the propagation of interference in time; the spatial causal edges connect the interference area nodes adjacent in space at the same time point, indicating the diffusion of interference in space; the source-area edges connect the interference source nodes and the interference area nodes, indicating the direct impact of human activities on specific areas. Finally, a directed graph containing all nodes and edges is constructed, providing a structural basis for causal relationship reasoning. The reason for the technical feature is that by clearly defining the nodes and edges, the causal network can intuitively represent the spatio-temporal propagation and causal sources of interference. The causal network structure provides a visual and systematic framework for causal reasoning, helping to identify complex causal chains.

[0120] Causal relationship reasoning:

[0121] After the causal network is constructed, causal relationship reasoning is carried out to evaluate and filter the edges in the causal network. Causal relationship 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 degree of the interference region nodes at adjacent time points. The calculation method is as follows: First, calculate the covariance of the interference flags at two time points, then calculate the standard deviations of the interference flags at the two time points respectively, and then divide the covariance by the product of the two standard deviations to obtain the correlation coefficient, which measures the continuity of the interference in time. If the calculated correlation coefficient is greater than the preset temporal threshold, the corresponding temporal causal edge is retained.

[0123] Spatial causal reasoning: The Moran's I index is used to evaluate the spatial correlation degree between adjacent grids at the same time point. The calculation method is as follows: First, calculate the deviation of the interference flags between the target grid and the neighborhood grids, multiply these deviations pairwise and sum them, and then normalize by the total deviation of the interference flags of all grids in the study area to obtain the Moran's I index, which measures the degree of aggregation or diffusion of the interference in space. 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: The Granger causality test is used to evaluate the causal influence of the interference source node on the interference region node. The calculation method is as follows: First, construct a time series model including the interference source activity and a time series model without the interference source activity, compare the goodness of fit of the two models, and calculate the p-value at the significance level. If the p-value is less than the preset significance level, the corresponding source-region edge is retained.

[0125] Using mature statistical methods such as the 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, providing a reliable basis for screening the edges in the causal network.

[0126] S5.3, Causal path analysis:

[0127] After completing the causal relationship reasoning, causal path analysis is performed to identify the influence paths of interference sources on the interference regions. The processing method is as follows: in the constructed and filtered causal network, for each interference region node, calculate the shortest path from all interference source nodes to this interference region node. The path weight is defined as the reciprocal of the association strength. The calculation process is as follows: starting from the interference source node, traverse step by step along the directed edges, 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 that has the greatest impact on this interference region node and its propagation path. Finally, count the causal sources of all interference region nodes, and aggregate to identify the main interference sources and their influence scopes.

[0128] The shortest path algorithm can efficiently identify the optimal causal chain between the interference source and the interference region, and reveal 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 artificial interference region map, geospatial data, and historical activity information generated in step S4, and realizes the analysis of the causal sources and propagation paths of the artificial interference regions. The network clearly reveals the impact of human activities as interference sources on natural resources and their spatio-temporal propagation laws. By causal path analysis, the main interference sources and their influence scopes are identified, providing a targeted basis for protection and intervention measures in the dynamic monitoring of natural resources.

[0130] The above formulas are all dimensionless and take their numerical values for calculation. The formulas are obtained by collecting a large amount of data for software simulation to obtain a formula closest to the actual situation. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0131] It should be noted that the system of the present invention can be deployed on the device itself to achieve embedded applications, or can also run on a PC with a user interface or other terminals, so as to meet various hardware environments and usage requirements.

[0132] Only some exemplary embodiments of the present invention have been described by way of illustration above. Undoubtedly, for those of ordinary skill in the art, without departing from the spirit and scope of the present invention, the described embodiments can be modified in various different ways. 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 text, relational terms such as first and second are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the term "comprising", "including" or any other variant thereof is intended to cover non-exclusive inclusion, so that a process, method, article or device comprising a series of elements not only includes those elements, but also includes other elements not expressly listed, or also includes elements inherent in such process, method, article or device. Without further limitation, an element defined by the statement "comprising an..." does not exclude the presence of additional identical elements in the process, method, article or device comprising the element.

[0134] As described above, the above is only the specific implementation manner of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present application can easily think of changes or substitutions, and all should be covered within the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claimed rights.

Claims

1. An unmanned aerial vehicle remote sensing Gaussian splash change detection method for dynamic monitoring of natural resources, characterized in that Including the steps: S1. Obtain multi-temporal UAV remote sensing point cloud data and convert it into a spatio-temporally continuous observation density field; S2. Conduct a time-axis gradient analysis on the observation density field, separate the low-frequency gradient component of natural evolution and the high-frequency gradient component of human-induced mutations, and generate a change pattern classification map; S3. Based on the low-frequency gradient component of natural evolution, obtain the predicted density field under the non-interference state; S4. Calculate the difference between the observation density field and the predicted density field at different spatial scales, perform consistency weighted fusion of the differences between scales into a comprehensive difference, dynamically set a judgment threshold in combination with the regional historical change amplitude and spatial heterogeneity, and locate the human interference area according to the comparison between the comprehensive difference and the dynamic threshold; S5. Combine geospatial data and historical activity information to construct a causal inference network, and analyze the causal sources and propagation paths of the human interference area.

2. The method for detecting Gaussian splash changes in UAV remote sensing for dynamic monitoring of natural resources according to claim 1, wherein Step S1 includes the following: Denoise the multi-temporal UAV remote sensing point cloud data, use a density-based clustering method to remove outliers to generate denoised point cloud data, align all denoised point cloud data to a unified geographic coordinate system through a registration method, divide the study area into uniform two-dimensional grids and subdivide them into multiple height layers in the vertical direction, allocate the denoised point cloud data to the corresponding two-dimensional grid cells and height layers according to the three-dimensional coordinates, calculate the local observation density value for the point set in each height layer using the Gaussian kernel density estimation method and generate a comprehensive observation density value through weighted summation, perform linear interpolation in the time dimension for the comprehensive observation density value to fill in the missing time-phase data and apply Gaussian filtering for spatio-temporal smoothing processing, and finally generate the observation density field.

3. The method for detecting Gaussian splash changes in UAV remote sensing for dynamic monitoring of natural resources according to claim 2, characterized in that, Step S2 includes the following: Conduct a time-axis gradient analysis on the observation density field data, calculate the density change gradient for the observation density value sequence of each spatial grid cell using the central difference method, perform a discrete Fourier transform on the density change gradient sequence to convert it into frequency-domain spectrum data, separate the spectrum data into low-frequency and high-frequency components according to a pre-set frequency threshold, perform an inverse discrete Fourier transform on the low-frequency and high-frequency components respectively to reconstruct the low-frequency and high-frequency gradient components in the time domain, calculate the relative intensity ratio of the high-frequency gradient component to the low-frequency gradient component and compare it with the classification threshold to generate the change pattern classification map.

4. The method for detecting Gaussian splash changes in UAV remote sensing for dynamic monitoring of natural resources according to claim 3, characterized in that, Step S3 includes the following: Based on the separated low-frequency gradient component, use the observation density value at the initial time point provided by the observation density field, reconstruct the predicted density value at subsequent time points by recursively accumulating the product of the low-frequency gradient component and the time interval, perform Gaussian filtering for smoothing processing on the reconstructed predicted density field in the time dimension, generate a smoothed predicted density value by performing Gaussian weighted averaging on the predicted density values of adjacent time points within the time window, and integrate the smoothed predicted density values into a time-space three-dimensional tensor to form the predicted density field under the non-interference state.

5. The method for detecting Gaussian splash changes in drone remote sensing for dynamic monitoring of natural resources according to claim 4, characterized in that, Step S4 includes the following: 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 as the difference value through Gaussian smoothing processing at multiple spatial scales. The consistency weighted fusion method is used to calculate the entropy value of the difference field in the local neighborhood at each spatial scale, and the weights are determined according to the entropy value, and then weighted average is performed to generate the comprehensive difference value. Combining the average absolute difference of the historical observed density value change of each spatial grid cell and the spatial heterogeneity calculated by the local Moran's I index, a dynamic threshold is set, and an artificial interference area map is generated by comparing the comprehensive difference value with the dynamic threshold.

6. The method for detecting Gaussian splash changes in drone remote sensing for dynamic monitoring of natural resources according to claim 5, wherein Step S5 includes the following: Integrate the artificial interference area map, geospatial data, and historical activity information and align them in the spatio-temporal dimension. Define the interference source nodes based on the historical activity information and the interference area nodes based on the artificial interference area map, and then construct a causal network including time causal edges, spatial causal edges, and source-region edges. Based on the filtered causal network, use the shortest path algorithm to identify the interference sources and their propagation paths, and count the main interference sources and the influence ranges.

7. The method for detecting Gaussian splash changes in drone remote sensing for dynamic monitoring of natural resources according to claim 6, characterized in that Step S5 also includes the following: Calculate the temporal correlation degree of the interference area nodes at adjacent time points through the Pearson correlation coefficient and compare it with the temporal threshold to filter the time causal edges.

8. The method for detecting Gaussian splash changes in UAV remote sensing for dynamic monitoring of natural resources according to claim 6, characterized in that, Step S5 also includes the following: Calculate the spatial correlation degree between adjacent grids at the same time point through the Moran's I index and compare it with the spatial threshold to filter the spatial causal edges.

9. The method for detecting Gaussian splash changes in UAV remote sensing for dynamic monitoring of natural resources according to claim 6, wherein Step S5 also includes the following: Calculate the causal influence of the interference source nodes on the interference area nodes through the Granger causality test and compare it with the significance level to filter the source-region edges.

Images