Intelligent analysis method for spatial-temporal variation of stoichiometric characteristics of forest plants

Abstract

The invention relates to the technical field of ecological environment monitoring and forest remote sensing, in particular to a forest plant stoichiometric characteristic spatial-temporal change intelligent analysis method, which comprises the following steps: S1, acquiring spectral reflectivity and element concentration values of plant leaves in a target area; s2, performing space-time registration processing; s3, constructing a three-dimensional space-time matrix, and generating a space-time cube with a dynamic weight identifier; s4, calculating the space-time gradient variation of the carbon-nitrogen-phosphorus ratio in each space unit; s5, performing feature decoupling on the space-time gradient variation; and S6, generating a stoichiometric characteristic three-dimensional atlas with altitude-time double coordinate axes. According to the method, by constructing the multi-factor driven three-dimensional space-time analysis model, decoupling analysis of long-term trend and short-term fluctuation of forest plant stoichiometric characteristics is achieved, the element flux mutation threshold value at the terrain turning position is accurately recognized, and the dynamic monitoring capacity of the ecological process is effectively improved.

Description

Technical Field

[0001] The present invention relates to the field of ecological environment monitoring and forest remote sensing technology, and in particular to an intelligent analysis method for the spatiotemporal changes of forest plant stoichiometric characteristics. Background Art

[0002] With the deepening of forest ecosystem research, the need for refined monitoring and analysis of the stoichiometric characteristics of elements such as carbon, nitrogen, and phosphorus in forest plant leaves has become increasingly urgent; hyperspectral remote sensing, ion mobility spectrometers, and multi-sensor fusion observation methods have been promoted in field surveys, which can provide a more comprehensive quantitative basis for the nutrient cycle, health status, and community succession of forest ecosystems; however, affected by the complexity of the terrain and the coupling of ecological factors, the observation data often have uneven sampling, temperature drift, environmental noise, and multi-dimensional interpolation problems at the spatiotemporal level, resulting in a lack of reliable and systematic analysis methods for the dynamic changes in the stoichiometric characteristics of plants.

[0003] Existing technologies have technical gaps in multi-source data fusion and correction of day-night temperature difference interference, making it difficult to accurately analyze carbon-nitrogen-phosphorus ratios under multiple environmental gradients. In addition, the interpolation methods for multidimensional environmental factors such as terrain slope, soil moisture content, and solar radiation are often single, and cannot take into account the separation of long-term trends and short-term fluctuations. In the absence of complete spatiotemporal registration and feature decoupling, researchers find it difficult to identify significant change points in element migration rates in forest ecosystems, and are unable to quantify element flux mutation thresholds at terrain turning points. In order to overcome the above shortcomings and meet the needs of accurate monitoring and early warning of forest ecological processes, an intelligent analysis method for spatiotemporal changes in forest plant stoichiometric characteristics is urgently needed to solve the problem. Summary of the invention

[0004] Based on the above objectives, the present invention provides an intelligent analysis method for the spatiotemporal changes of forest plant stoichiometric characteristics.

[0005] The intelligent analysis method of spatiotemporal changes of forest plant stoichiometric characteristics includes the following steps: S1: Use a multispectral drone equipped with an ion mobility spectrometer to synchronously obtain the spectral reflectance, element concentration values, and corresponding longitude and latitude coordinates, altitude, and acquisition timestamp of plant leaves in the target area to form an original multi-source data set; S2: Perform spatiotemporal registration processing on the original multi-source data set, use the improved Lomb-Scargle periodic algorithm to eliminate the signal drift caused by the day and night temperature difference, and output the standardized multi-source data; S3: Construct a three-dimensional space-time matrix, perform spatial interpolation and reorganization of the standardized multi-source data according to the three dimensions of elevation gradient, solar radiation intensity, and soil moisture content, and generate a space-time cube with dynamic weight labels; S4: Based on the dynamic weight identification in the space-time cube, a nonlinear mapping relationship between terrain slope and element migration rate is established, and the spatiotemporal gradient change of carbon-nitrogen-phosphorus ratio in each spatial unit is calculated through an adaptive sliding window; S5: Decouple the temporal and spatial gradient changes to separate the long-term trend component dominated by geographical factors and the short-term fluctuation component caused by environmental mutations; S6: Based on the decoupled trend component and fluctuation component, a three-dimensional map of stoichiometric characteristics with dual coordinate axes of altitude and time is generated, and the element flux mutation threshold at the terrain turning point is marked.

[0006] Optionally, the S1 specifically includes: S11: A multispectral camera array and an ion mobility spectrometer are integrated in the belly of the drone. The multispectral camera includes optical sensors with five wavelengths of 450nm, 550nm, 670nm, 780nm and 850nm. The inlet of the ion mobility spectrometer is equipped with a micro air pump to maintain a 0.5mm distance from the surface of the plant leaf. S12: When the UAV hovers at a height of 3m above the target plant, the multispectral camera and the ion mobility spectrometer are triggered by the laser rangefinder to start acquisition synchronously. The laser rangefinder sends a trigger pulse signal when the error of the distance to the top of the canopy is less than ±0.2m. S13: The RTK-GPS module on the drone records the longitude and latitude coordinates in real time. The barometric altimeter updates the altitude value every 0.5 seconds. The Beidou timing chip marks the acquisition timestamp with 50ms accuracy. S14: The multispectral camera completes the spectral reflectance acquisition of five bands within 20ms after being triggered. At the same time, the ion mobility spectrometer obtains the concentration values ​​of the three elements C, N, and P within a 150ms detection cycle. The two are synchronized at the hardware level through the FPGA chip. S15: Encapsulate the spectral reflectance data, element concentration values, geographic location data and timestamp according to a unified time reference to form an original multi-source data set.

[0007] Optionally, the S2 specifically includes: S21: sort the spectral reflectance and element concentration values ​​in the original multi-source data set formed in S1 in time series, extract the single plant observation sequence at each time point and unify it into a sampling interval of one group every 10 minutes; S22: Using the Lomb-Scargle periodic analysis method, we construct a power spectrum distribution diagram in the frequency domain for data sampled at unequal time intervals, and set 0.0417 Hz as the main periodic reference frequency to identify the periodic signal fluctuation components caused by the day-night temperature difference; S23: The sampling weight factor and temperature covariate are introduced to improve the standard Lomb-Scargle algorithm and obtain the modified weighted power spectrum; S24: Based on the weighted periodic model, the temperature-driven signal drift sequence is reconstructed, and the drift trend curve of each observation point is fitted using the least squares method, and then the trend is removed point by point from the original spectral reflectance and element concentration values ​​to obtain standardized data after temperature drift correction; S25: The drift-corrected spectral reflectance and element concentration values ​​are re-bound to the latitude and longitude coordinates, altitude and timestamp to form a standardized multi-source data set.

[0008] Optionally, the S24 specifically includes: S241: Based on the weighted power spectrum of S23, extract the amplitude coefficients of the cosine term and the sine term at the reference frequency of the corresponding main cycle, and reconstruct the drift signal sequence of each observation point. The reconstruction expression is: ,in, For the Drift estimate for each observation point; is the drift signal amplitude; is the reference angular frequency; For the observation time; is the drift phase angle; S242: Use the least squares method to fit the drift sequence and construct the error square function , whose expression is: ,in, For the observations, representing raw spectral reflectance or element concentration values; S243: Perform signal correction operation on each observation point, remove the estimated drift trend from the original observation value, and obtain the drift-corrected standardized data, which is expressed as: ,in, are normalized values ​​after correction for temperature driven drift.

[0009] Optionally, the S3 specifically includes: S31: Based on the standardized multi-source data obtained in S2, the altitude data corresponding to each observation data point is , construct the initial data sequence along the elevation direction, and perform linear interpolation on the elevation sequence to make the spatial interval of the data sequence on the elevation axis uniform and fixed. ; S32: Find the corresponding satellite remote sensing solar radiation intensity data according to the latitude and longitude coordinates and timestamp of the location of the observation data point , calculate the solar radiation intensity value of each data point, and use the inverse distance weighted interpolation method to spatially resample the solar radiation intensity to generate a two-dimensional spatial continuous distribution map of solar radiation intensity; S33: Use soil moisture sensor to collect soil moisture value in target area , and perform spatial Kriging interpolation calculation based on the latitude and longitude positions of the soil moisture sampling points to generate a spatial continuous distribution map of soil moisture; S34: Based on the data sequence after interpolation in the elevation direction, the solar radiation intensity data obtained by S32 interpolation and the soil moisture data obtained by S33 interpolation are superimposed to form a three-dimensional spatial data structure, and a dynamic weight identification value is assigned to each spatial grid unit. ; S35: The three-dimensional spatial data structure reconstructed in S34 and its dynamic weight identifier are uniformly bound to the corresponding latitude and longitude coordinates, elevation position and timestamp to form a space-time cube with a dynamic weight identifier.

[0010] Optionally, the S4 specifically includes: S41: Based on the space-time cube formed in S3, the terrain slope value of each spatial unit is calculated according to the latitude and longitude coordinates and elevation values ​​of the spatial unit, and the initial mapping relationship between the slope and the dynamic weight identifier is determined in combination with the dynamic weight identifier value corresponding to the spatial unit; S42: Calculate the migration rate of the elements in each spatial unit according to the change sequence of the element concentration value in the spatial unit over time, and weight the migration rate using the dynamic weight identifier to establish a mapping relationship between the migration rate and the dynamic weight identifier; S43: Based on the initial mapping relationship between slope and weight and the mapping relationship between migration rate and weight, a nonlinear mapping function between terrain slope and element migration rate is established by using a nonlinear regression method; S44: construct an adaptive sliding window with each space unit in the space-time cube as the center, and the size of the sliding window is adaptively adjusted according to the dynamic weight identifier of the central unit. The larger the weight identifier value, the smaller the window size; S45: within the adaptive sliding window, respectively counting the average concentrations of carbon, nitrogen and phosphorus, and calculating the carbon, nitrogen and phosphorus concentration ratio corresponding to the central unit of the window; S46: Slide the adaptive sliding window gradually along the time axis and the space axis, differentiate the carbon-nitrogen-phosphorus ratios of adjacent moments and adjacent spatial units, calculate the spatiotemporal gradient change of the carbon-nitrogen-phosphorus ratio corresponding to each spatial unit, and generate a spatiotemporal gradient change data sequence of the carbon-nitrogen-phosphorus ratio.

[0011] Optionally, the expression of the nonlinear mapping function in S43 is: ,in, is the predicted value of element migration rate; is the terrain slope value; is the saturation value parameter of the element migration rate, indicating the theoretical maximum value of the element migration rate; is the nonlinear change rate parameter of the curve; is the turning point parameter of the slope response, indicating the slope value corresponding to the transition from low speed to high speed of element migration rate; is the natural logarithm base constant.

[0012] Optionally, the S5 specifically includes: S51: sorting the spatiotemporal gradient change data sequence of the carbon-nitrogen-phosphorus ratio calculated in S4, extracting the gradient change value of each spatial unit in the continuous time series, and forming a time series feature sequence; S52: The Daubechies wavelet basis function is used to perform three-layer discrete wavelet decomposition on the gradient change time series of each spatial unit, decomposing the original sequence into low-frequency approximate components and high-frequency detail components, where the low-frequency component represents long-term trend information and the high-frequency component represents mutation fluctuation information; S53: retain the low-frequency approximate component, remove the noise interference term in the high-frequency detail component, and restore the long-term trend characteristic curve through inverse wavelet transform to obtain the long-term trend component corresponding to each spatial unit; S54: for the high-frequency detail component, performing correlation analysis between spatially adjacent units, and constructing a distribution map of the high-frequency component in three-dimensional space; S55: Using spatial coordinates as variables, spatial Kriging interpolation is performed on the high-frequency detail components to estimate the propagation path and intensity of sudden fluctuations in space, and to construct a spatially continuous short-term fluctuation component field.

[0013] Optionally, the S55 specifically includes: S551: Based on the high-frequency detail component distribution map of the spatial unit obtained in S54, the spatial unit with a significant fluctuation peak is taken as the initial source point to determine the starting position of the spatial propagation of the mutation fluctuation; S552: Taking the source point as the spatial interpolation center, define the interpolation influence radius, extract the high-frequency detail component values ​​of adjacent spatial units within the radius, and use the variation function method to calculate the structural correlation degree of the high-frequency detail component in space, and determine the dominant direction of spatial propagation of the high-frequency detail component; S553: ​​Perform spatial Kriging interpolation step by step along the determined spatial propagation advantage direction to obtain a continuous prediction field of the spatial distribution of high-frequency detail components, and identify the fluctuation propagation path according to the spatial variation trend of the prediction field value; S554: Calculate the wave propagation intensity based on the spatial interpolation result of the high-frequency detail component, determine the numerical difference of the high-frequency detail component between adjacent spatial units on the propagation path, calculate the fluctuation intensity change gradient, and clearly indicate the degree of change of the fluctuation intensity on the propagation path; S555: Integrate the wave propagation path and the wave propagation intensity gradient change to form a complete mutation wave space propagation path map and intensity change field, and output the spatially continuous short-term wave component.

[0014] Optionally, the S6 specifically includes: S61: Based on the long-term trend component obtained by decoupling in S5, the initial three-dimensional surface of the long-term trend characteristics is constructed by using the three-dimensional visualization method, with the time axis and altitude axis as the horizontal and vertical coordinates and the long-term trend component value as the vertical coordinate; S62: Using the short-term fluctuation component obtained in S5, with the same three-dimensional coordinate axis system, by superimposing the short-term fluctuation values, a complete three-dimensional map of the trend component and the fluctuation component is formed; S63: Perform spatial derivative analysis on the three-dimensional map after the trend component and the fluctuation component are superimposed, calculate the derivative value of the element flux in the map with the change of altitude, determine the spatial position where the derivative in the altitude direction has a significant mutation, and mark it as the terrain turning point; S64: Taking the determined terrain turning point as the center, analyzing the element flux change trend at the location along the time axis, calculating the flux change rate, and setting a significance threshold, thereby determining the element flux mutation threshold; S65: The determined terrain turning points and the corresponding element flux mutation thresholds are clearly marked on the three-dimensional map and visualized through obvious graphic labels.

[0015] Beneficial effects of the present invention: The present invention, by fusing multi-sensor observation data with an improved diurnal temperature difference drift correction algorithm, constructs a three-dimensional space-time cube and adopts a composite algorithm combining wavelet transform and spatial Kriging interpolation, effectively integrating the influence of multi-dimensional ecological factors on the stoichiometric characteristics of forest plants; the scheme significantly improves the accuracy and resolution of space-time data under multi-factor gradients such as elevation, solar radiation and soil moisture content, enables the simultaneous identification of long-term trends and short-term fluctuations, and provides a clear space-time analysis capability for the dynamic evolution of the carbon-nitrogen-phosphorus ratio; The present invention uses a nonlinear mapping model to accurately characterize the quantitative relationship between terrain slope and element migration rate, and realizes the significance determination of element flux through derivative analysis and threshold calculation. This scheme greatly improves the accuracy and reliability of the stoichiometric characteristic analysis of forest ecosystems. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] In order to more clearly illustrate the technical solutions in the present invention or the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings in the following description are only for the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0017] Figure 1 A schematic diagram of an intelligent analysis method according to an embodiment of the present invention; Figure 2 Schematic diagram of the calculation process of the spatiotemporal gradient variation according to an embodiment of the present invention. DETAILED DESCRIPTION

[0018] The present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. At the same time, it is explained here that in order to make the embodiments more detailed, the following embodiments are the best and preferred embodiments, and those skilled in the art may also adopt other alternatives to implement some known technologies; and the accompanying drawings are only for more specific description of the embodiments, and are not intended to specifically limit the present invention.

[0019] It should be noted that the references to "one embodiment", "embodiment", "exemplary embodiments", "some embodiments" and the like in the specification indicate that the embodiments described may include specific features, structures or characteristics, but not every embodiment may include the specific features, structures or characteristics. In addition, when a specific feature, structure or characteristic is described in conjunction with an embodiment, it should be within the knowledge of a person skilled in the art to implement such feature, structure or characteristic in conjunction with other embodiments (whether or not explicitly described).

[0020] In general, a term can be understood, at least in part, from its use in context. For example, depending, at least in part, on the context, the term "one or more" as used herein can be used to describe any feature, structure, or characteristic in the singular sense, or can be used to describe a combination of features, structures, or characteristics in the plural sense. Additionally, the term "based on" can be understood as not necessarily intended to convey an exclusive set of factors, but can instead, depending, at least in part, on the context, allow for the presence of other factors that are not necessarily explicitly described.

[0021] like Figure 1-Figure 2 As shown, the intelligent analysis method of spatiotemporal changes of forest plant stoichiometric characteristics includes the following steps: S1: Use a multispectral drone equipped with an ion mobility spectrometer to synchronously obtain the spectral reflectance, element concentration values, and corresponding longitude and latitude coordinates, altitude, and acquisition timestamp of plant leaves in the target area to form an original multi-source data set; S2: Perform spatiotemporal registration processing on the original multi-source data set, use the improved Lomb-Scargle periodic algorithm to eliminate the signal drift caused by the day and night temperature difference, and output the standardized multi-source data; S3: Construct a three-dimensional space-time matrix, perform spatial interpolation and reorganization of the standardized multi-source data according to the three dimensions of elevation gradient, solar radiation intensity, and soil moisture content, and generate a space-time cube with dynamic weight labels; S4: Based on the dynamic weight identification in the space-time cube, a nonlinear mapping relationship between terrain slope and element migration rate is established, and the spatiotemporal gradient change of carbon-nitrogen-phosphorus ratio in each spatial unit is calculated through an adaptive sliding window; S5: Decouple the temporal and spatial gradient changes to separate the long-term trend component dominated by geographical factors and the short-term fluctuation component caused by environmental mutations; S6: Based on the decoupled trend component and fluctuation component, a three-dimensional map of stoichiometric characteristics with dual coordinate axes of altitude and time is generated, and the element flux mutation threshold at the terrain turning point is marked.

[0022] S1 specifically includes: S11: A multispectral camera array and an ion mobility spectrometer are integrated in the belly of the drone. The multispectral camera includes optical sensors with five wavelengths of 450nm, 550nm, 670nm, 780nm and 850nm. The inlet of the ion mobility spectrometer is equipped with a micro air pump to maintain a 0.5mm distance from the surface of the plant leaf. S12: When the UAV hovers at a height of 3m above the target plant, the multispectral camera and the ion mobility spectrometer are triggered by the laser rangefinder to start acquisition synchronously. The laser rangefinder sends a trigger pulse signal when the error of the distance to the top of the canopy is less than ±0.2m. S13: The RTK-GPS module on the drone records the longitude and latitude coordinates in real time. The barometric altimeter updates the altitude value every 0.5 seconds. The Beidou timing chip marks the acquisition timestamp with 50ms accuracy. S14: The multispectral camera completes the spectral reflectance acquisition of five bands within 20ms after being triggered. At the same time, the ion mobility spectrometer obtains the concentration values ​​of the three elements C, N, and P within a 150ms detection cycle. The two are synchronized at the hardware level through the FPGA chip. S15: Encapsulate the spectral reflectance data, element concentration values, geographic location data and timestamp according to a unified time reference to form an original multi-source data set; the above steps integrate a multispectral camera and an ion mobility spectrometer on the UAV platform, and introduce laser ranging, RTK-GPS, barometric altimeter and BeiDou timing technology to achieve high-temporal and spatial resolution in-situ synchronous acquisition of the chemical metrological characteristics of forest plant leaves, providing a high-precision, multi-dimensional data foundation for subsequent spatiotemporal analysis.

[0023] S2 specifically includes: S21: sort the spectral reflectance and element concentration values ​​in the original multi-source data set formed in S1 in time series, extract the single plant observation sequence at each time point and unify it into a sampling interval of one group every 10 minutes; S22: Using the Lomb-Scargle periodic analysis method, we construct a power spectrum distribution diagram in the frequency domain for data sampled at different time intervals, and set 0.0417 Hz (i.e., 24-hour period) as the main periodic reference frequency to identify the periodic signal fluctuation components caused by the day-night temperature difference; S23: Introducing sampling weight factors and temperature covariates to improve the standard Lomb-Scargle algorithm and obtain the corrected weighted power spectrum , where the weight factor is calculated from the ambient temperature change rate, and the covariate is the temperature curve obtained by the meteorological station in the same period to improve the sensitivity and accuracy of period analysis; The standard Lomb-Scargle algorithm is improved as follows: S231: Set the temperature change rate covariate to , whose expression is: ,in, For the moment External ambient temperature; S232: Based on Calculate the sampling weight factor for each observation point , defined as follows: ,in, For the The sampling time, is the sensitivity control coefficient, is the absolute value of the temperature change rate at that moment; S233: Introducing the Lomb-Scargle algorithm, the original observations Weighted, modified weighted power spectrum The calculation formula is: ,in, For the An observation value represents a single measured spectral reflectance or element concentration; is the weighted mean, defined as ; is the angular frequency; is the phase offset correction amount.

[0024] S24: Based on the weighted periodic model, the temperature-driven signal drift sequence is reconstructed, and the drift trend curve of each observation point is fitted using the least squares method, and then the trend is removed point by point from the original spectral reflectance and element concentration values ​​to obtain standardized data after temperature drift correction; S25: The drift-corrected spectral reflectance and element concentration values ​​are re-bound to the longitude and latitude coordinates, altitude and timestamp to form a standardized multi-source data set as the input basis for subsequent spatiotemporal analysis steps; the above steps introduce temperature weights and covariate correction mechanisms to the Lomb-Scargle periodic algorithm, combined with reverse elimination of drift trends, to improve the algorithm's modeling accuracy for diurnal temperature-driven errors, ensure the comparability and stability of multi-source observation data in the time dimension, and lay a solid data foundation for the subsequent construction of spatial interpolation structure and feature decoupling.

[0025] S24 specifically includes: S241: Based on the weighted power spectrum of S23, extract the amplitude coefficients of the cosine term and the sine term at the reference frequency of the corresponding main cycle, and reconstruct the drift signal sequence of each observation point. The reconstruction expression is: ,in, For the Drift estimate for each observation point; is the drift signal amplitude; is the reference angular frequency, corresponding to the diurnal cycle; For the observation time; is the drift phase angle; S242: Use the least squares method to fit the drift sequence and construct the error square function , whose expression is: ,in, For the observations, representing raw spectral reflectance or element concentration values; S243: Perform signal correction operation on each observation point, remove the estimated drift trend from the original observation value, and obtain the drift-corrected standardized data, which is expressed as: ,in, It is the standardized value after temperature-driven drift correction. By constructing a weighted periodic model and combining it with the least squares fitting method, the drift component driven by temperature is accurately reconstructed and the trend is eliminated point by point, which effectively improves the stability and physical consistency of the observation data, so that subsequent spatiotemporal interpolation and multi-factor modeling have a unified baseline reference condition.

[0026] S3 specifically includes: S31: Based on the standardized multi-source data obtained in S2, the altitude data corresponding to each observation data point is , construct the initial data sequence along the elevation direction, and perform linear interpolation on the elevation sequence to make the spatial interval of the data sequence on the elevation axis uniform and fixed. ; S32: Find the corresponding satellite remote sensing solar radiation intensity data according to the latitude and longitude coordinates and timestamp of the location of the observation data point , calculate the solar radiation intensity value of each data point, and use the inverse distance weighted interpolation method to spatially resample the solar radiation intensity to generate a two-dimensional spatial continuous distribution map of solar radiation intensity; The calculation formula for inverse distance weighted interpolation is: ,in, For spatial location The interpolated solar radiation intensity value at ; For the The solar radiation intensity values ​​at adjacent observation data points; To be interpolated position to The Euclidean distance of each data point in meters; : Distance attenuation coefficient, the value range is ; The number of adjacent data points used for interpolation calculation is fixed at 8; S33: Use soil moisture sensor to collect soil moisture value in target area , and perform spatial Kriging interpolation calculation based on the latitude and longitude positions of the soil moisture sampling points to generate a spatial continuous distribution map of soil moisture; The Kriging interpolation formula is: ,in: For location The soil moisture content after interpolation at For the The soil moisture content measured at each sampling point; is the weight coefficient determined by the Kriging algorithm optimization; The number of adjacent sampling points used in Kriging interpolation calculation is fixed at 12; S34: Based on the data sequence after interpolation in the elevation direction, the solar radiation intensity data obtained by S32 interpolation and the soil moisture data obtained by S33 interpolation are superimposed to form a three-dimensional spatial data structure, and a dynamic weight identification value is assigned to each spatial grid unit. , whose expression is: ,in, It is a dimensionless index with a dynamic weight identification value ranging from 0 to 1; is the altitude value of the spatial unit; is the solar radiation intensity value of the space unit; is the soil moisture value of the spatial unit; is the maximum elevation in the study area; is the maximum solar radiation intensity value in the study area; is the maximum soil moisture value in the study area; S35: The three-dimensional spatial data structure reconstructed in S34 and its dynamic weight identifier are uniformly bound to the corresponding longitude and latitude coordinates, elevation position and timestamp to form a space-time cube with a dynamic weight identifier; through the above interpolation and reconstruction steps, a clear multi-factor spatial interpolation method is used to effectively integrate multidimensional environmental factors such as elevation gradient, solar radiation intensity, and soil moisture content, so that the reconstructed three-dimensional data structure has clear spatiotemporal continuity and high environmental representativeness. At the same time, the introduction of dynamic weight indicators ensures that different spatial units have clear differences in ecological characteristic analysis, providing a rigorous and scientific data basis for subsequent nonlinear relationship modeling and trend decoupling analysis.

[0027] S4 specifically includes: S41: Based on the space-time cube formed in S3, the terrain slope value of each spatial unit is calculated according to the latitude and longitude coordinates and elevation values ​​of the spatial unit, and the initial mapping relationship between the slope and the dynamic weight identifier is determined in combination with the dynamic weight identifier value corresponding to the spatial unit; The terrain slope value calculation formula is: ,in: is the terrain slope value, in radians; is the elevation difference between the target space unit and the adjacent space unit, in meters; is the horizontal distance between the target space unit and the adjacent space unit, in meters; The method for determining the initial slope-weight mapping relationship in S41 is as follows: Firstly, data pairs are constructed based on the slope values ​​of all spatial units in the space-time cube and the corresponding dynamic weight identification values ​​to form a training sample set; Then, based on the spline function, determine the slope value Identify values ​​for independent variables and dynamic weights is the initial nonlinear mapping relationship of the dependent variable. The cubic spline interpolation method is used in the specific calculation, in which the number and position of nodes are determined by the uniform distribution of data point density; Finally, the least squares criterion is used to optimize the parameters in the spline function so that the sum of squares of the errors between the predicted dynamic weight identification value and the actual observed value is minimized.

[0028] S42: Calculate the migration rate of the elements in each spatial unit according to the change sequence of the element concentration value in the spatial unit over time, and weight the migration rate using the dynamic weight identifier to establish a mapping relationship between the migration rate and the dynamic weight identifier; The element migration rate calculation formula is: ,in, is the element migration rate, the weighted concentration change rate, in milligrams per kilogram per hour; For an element at time The concentration value is in milligrams per kilogram; for the same element at the previous moment The concentration value is in milligrams per kilogram; is the time interval between two consecutive observation times, in hours; is the dynamic weight identification value of the target spatial unit; The method for determining the migration rate-weight mapping relationship in the above S42 is as follows: First, calculate the element migration rate data obtained from multiple observations in each spatial unit and the dynamic weight identification value of the corresponding unit Constructing data pairs to generate training data sets of element migration rates and dynamic weight identifiers; Then, identify the value with dynamic weight is the independent variable, element migration rate As the dependent variable, an exponential function is used for fitting, and the corresponding fitting parameters in the exponential function are optimized using the nonlinear least square method to minimize the sum of squares of the residuals between the element migration rate predicted by the function and the actual migration rate observations. Finally, based on the fitted exponential function model, the quantitative impact of changes in dynamic weight identification values ​​on element migration rates can be clearly described.

[0029] S43: Based on the initial mapping relationship between slope and weight and the mapping relationship between migration rate and weight, a nonlinear mapping function between terrain slope and element migration rate is established by using a nonlinear regression method; S44: construct an adaptive sliding window with each space unit in the space-time cube as the center, and the size of the sliding window is adaptively adjusted according to the dynamic weight identifier of the central unit. The larger the weight identifier value, the smaller the window size; S45: within the adaptive sliding window, respectively counting the average concentrations of carbon, nitrogen and phosphorus, and calculating the carbon, nitrogen and phosphorus concentration ratio corresponding to the central unit of the window; S46: Slide the adaptive sliding window gradually along the time axis and the space axis, differentiate the carbon-nitrogen-phosphorus ratios of adjacent moments and adjacent spatial units, calculate the spatiotemporal gradient change of the carbon-nitrogen-phosphorus ratio corresponding to each spatial unit, and generate a spatiotemporal gradient change data sequence of the carbon-nitrogen-phosphorus ratio; the above steps establish an accurate nonlinear mapping relationship between terrain slope and element migration rate by utilizing dynamic weight markers, and further introduce an adaptive sliding window mechanism, which effectively solves the analysis bias caused by a single spatial scale, improves the ability to analyze the element migration process and the spatiotemporal gradient of the carbon-nitrogen-phosphorus ratio change, and lays a reliable data support for the precise decoupling analysis of subsequent trends and fluctuation characteristics.

[0030] The expression of the nonlinear mapping function in S43 is: ,in, is the predicted value of element migration rate, in milligrams per kilogram per hour; is the terrain slope value; It is the saturation value parameter of element migration rate, indicating the theoretical maximum value of element migration rate, in milligrams per kilogram per hour; is the nonlinear change rate parameter of the curve, which determines the sensitivity of the slope to the change of migration rate; is the turning point parameter of the slope response, indicating the slope value corresponding to the transition from low speed to high speed of the element migration rate, in radians; It is the natural logarithm base constant (approximately 2.71828), dimensionless.

[0031] S5 specifically includes: S51: sorting the spatiotemporal gradient change data sequence of the carbon-nitrogen-phosphorus ratio calculated in S4, extracting the gradient change value of each spatial unit in the continuous time series, and forming a time series feature sequence; S52: The Daubechies wavelet basis function is used to perform three-layer discrete wavelet decomposition on the gradient change time series of each spatial unit, decomposing the original sequence into low-frequency approximate components and high-frequency detail components, where the low-frequency component represents long-term trend information and the high-frequency component represents mutation fluctuation information; The calculation formula of three-layer discrete wavelet decomposition is: ,in, It is the time series of the gradient change of the original carbon-chlorine-phosphorus ratio; It is the low-frequency approximate component obtained after the least-layer decomposition, representing the long-term trend information; For the The high-frequency detail components after layer decomposition represent the mutation fluctuation information of different scales; is the number of discrete wavelet decomposition layers, set to 3; is the time series index; S53: retain the low-frequency approximate component, remove the noise interference term in the high-frequency detail component, and restore the long-term trend characteristic curve through inverse wavelet transform to obtain the long-term trend component corresponding to each spatial unit; The inverse wavelet transform formula is: ,in, It is the long-term trend characteristic curve after recovery; For the Wavelet approximation coefficients at layer scale; For the Daubechies wavelet scaling function corresponding to the layer scale; is the scale coefficient index; S54: for the high-frequency detail component, performing correlation analysis between spatially adjacent units, and constructing a distribution map of the high-frequency component in three-dimensional space; The formula for spatial correlation analysis is: ,in, For space unit With adjacent space units The spatial correlation coefficient of the high-frequency detail components between them; Space Unit With space unit In the The value of the high-frequency detail component at each time point; Space Unit With space unit The temporal mean of the high-frequency detail component sequence; is the total number of time series data involved in the calculation; S55: Taking spatial coordinates as variables, spatial Kriging interpolation is performed on the high-frequency detail components to estimate the propagation path and intensity of mutation fluctuations in space, and to construct a spatially continuous short-term fluctuation component field; the above steps combine wavelet transform with spatial Kriging interpolation to perform dual-domain decoupling of the spatiotemporal gradient changes of the carbon-nitrogen-phosphorus ratio, and extract the trend terms dominated by geographical factors and the mutation terms driven by environmental disturbances respectively. This not only improves the accuracy of feature separation, but also enhances the operability of ecological anomaly identification and mechanism tracing, providing complete data support for terrain turning point anomaly identification and map visualization.

[0032] S55 specifically includes: S551: Based on the high-frequency detail component distribution map of the spatial unit obtained in S54, the spatial unit with a significant fluctuation peak is taken as the initial source point to determine the starting position of the spatial propagation of the mutation fluctuation; S552: Taking the source point as the spatial interpolation center, define the interpolation influence radius, extract the high-frequency detail component values ​​of adjacent spatial units within the radius, and use the variation function method to calculate the structural correlation degree of the high-frequency detail component in space, and determine the dominant direction of spatial propagation of the high-frequency detail component; The variogram method includes the following sub-steps: Taking the source point as the center, determine multiple spatial sampling point pairs in different directions; Calculate the square mean of the difference between the high-frequency detail components of each pair of spatial units in each direction; The direction with the smallest square mean is selected as the dominant direction of spatial propagation of high-frequency detail components.

[0033] S553: ​​Perform spatial Kriging interpolation step by step along the determined spatial propagation advantage direction to obtain a continuous prediction field of the spatial distribution of high-frequency detail components, and identify the fluctuation propagation path according to the spatial variation trend of the prediction field value; S554: Calculate the wave propagation intensity based on the spatial interpolation result of the high-frequency detail component, determine the numerical difference of the high-frequency detail component between adjacent spatial units on the propagation path, calculate the fluctuation intensity change gradient, and clearly indicate the degree of change of the fluctuation intensity on the propagation path; S555: The fluctuation propagation path and the fluctuation propagation intensity gradient change are integrated to form a complete mutation fluctuation spatial propagation path map and intensity change field, and output the spatially continuous short-term fluctuation component; through the above estimation steps, the continuous expression of mutation fluctuation information in the spatial dimension and the identification of intensity changes are accurately realized, so that the analysis of short-term mutation fluctuations has clear propagation trajectories and quantitative intensity information, which effectively enhances the ecological interpretability and predictive analysis accuracy of the mutation fluctuation process, and provides a reliable data basis for the subsequent accurate identification of terrain turning points and ecological early warning.

[0034] S6 specifically includes: S61: Based on the long-term trend component obtained by decoupling in S5, the initial three-dimensional surface of the long-term trend characteristics is constructed by using the three-dimensional visualization method, with the time axis and altitude axis as the horizontal and vertical coordinates and the long-term trend component value as the vertical coordinate; S62: Using the short-term fluctuation components obtained in S5, with the same three-dimensional coordinate axis system, by superimposing the short-term fluctuation values, a complete three-dimensional map of the trend component and the fluctuation component is formed, which is used to clearly reflect the comprehensive change trend and fluctuation of the element stoichiometric characteristics in the altitude and time dimensions; S63: Perform spatial derivative analysis on the three-dimensional map after the trend component and the fluctuation component are superimposed, calculate the derivative value of the element flux in the map with the change of altitude, determine the spatial position where the derivative in the altitude direction has a significant mutation, and mark it as the terrain turning point; The spatial derivative calculation formula is: ,in, is the derivative of the element flux in the altitude direction, in milligrams per kilogram per meter; Overlay values ​​for trends and fluctuations of element fluxes in three-dimensional maps; : Altitude, in meters; Time, in hours; S64: Taking the determined terrain turning point as the center, analyzing the element flux change trend at the location along the time axis, calculating the flux change rate, and setting a significance threshold, thereby determining the element flux mutation threshold; The calculation formula of element flux change rate is: ,in: is the element flux change rate, in milligrams per kilogram per hour; Indicates the corresponding altitude at the terrain turning point Place, time The element flux value of Indicates the same location , the previous moment The element flux value of is the time interval between two consecutive observation times, in hours; The steps to determine the element flux mutation threshold are as follows: S641: At the altitude corresponding to the terrain turning point, collect the element flux change data at multiple consecutive observation times, and calculate the mean and standard deviation of these data to obtain the basic statistical characteristics of the element flux fluctuation at this location; S642: Set the significance level α, and calculate the significance threshold that meets the α level based on the critical value of the normal distribution or other applicable distribution, which is used to measure whether the element flux change exceeds the normal fluctuation range; S643: Compare the element flux change rate at each moment at the altitude with the significance threshold. If the flux change rate at a certain moment exceeds the significance threshold, it is determined that a sudden change in the element flux occurs at that moment.

[0035] S65: The determined terrain turning points and the corresponding element flux mutation thresholds are clearly marked on the three-dimensional map, and visualized through obvious graphic symbols, so as to form a three-dimensional map of stoichiometric characteristics with altitude-time dual coordinate axes that is clearly marked and has significant ecological significance; the above steps use precise three-dimensional spatial derivative analysis and significance threshold determination methods to enable the generated three-dimensional map of stoichiometric characteristics to intuitively and clearly reveal the abnormal change points of element flux in the time and space dimensions, thereby strengthening the dynamic monitoring of nutrient migration characteristics of forest ecosystems and the early warning analysis capabilities of ecological processes.

[0036] The present invention covers any substitution, modification, equivalent method and scheme made on the essence and scope of the present invention. In order to make the public have a thorough understanding of the present invention, specific details are described in detail in the following preferred embodiments of the present invention, but those skilled in the art can fully understand the present invention without the description of these details. In addition, in order to avoid unnecessary confusion about the essence of the present invention, well-known methods, processes, procedures, components and circuits are not described in detail.

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

Claims

1. An intelligent analysis method for the spatiotemporal variation of forest plant stoichiometric characteristics, characterized in that: The following steps are involved: S1: Use a multispectral drone equipped with an ion mobility spectrometer to synchronously obtain the spectral reflectance, element concentration values, and corresponding longitude and latitude coordinates, altitude, and acquisition timestamp of plant leaves in the target area to form an original multi-source data set; S2: Perform spatiotemporal registration processing on the original multi-source data set, use the improved Lomb-Scargle periodic algorithm to eliminate the signal drift caused by the day and night temperature difference, and output the standardized multi-source data; S3: Construct a three-dimensional space-time matrix, perform spatial interpolation and reorganization of the standardized multi-source data according to the three dimensions of elevation gradient, solar radiation intensity, and soil moisture content, and generate a space-time cube with dynamic weight labels; S4: Based on the dynamic weight identification in the space-time cube, a nonlinear mapping relationship between terrain slope and element migration rate is established, and the spatiotemporal gradient change of carbon-nitrogen-phosphorus ratio in each spatial unit is calculated through an adaptive sliding window; S5: Decouple the temporal and spatial gradient changes to separate the long-term trend component dominated by geographical factors and the short-term fluctuation component caused by environmental mutations; S6: Based on the decoupled trend component and fluctuation component, a three-dimensional map of stoichiometric characteristics with dual coordinate axes of altitude and time is generated, and the element flux mutation threshold at the terrain turning point is marked.

2. The method for intelligent analysis of spatiotemporal changes in forest plant stoichiometric characteristics according to claim 1 is characterized in that: The S1 specifically includes: S11: A multispectral camera array and an ion mobility spectrometer are integrated in the belly of the drone. The multispectral camera includes optical sensors with five wavelengths of 450nm, 550nm, 670nm, 780nm and 850nm. The inlet of the ion mobility spectrometer is equipped with a micro air pump to maintain a 0.5mm distance from the surface of the plant leaf. S12: When the UAV hovers at a height of 3m above the target plant, the multispectral camera and the ion mobility spectrometer are triggered by the laser rangefinder to start acquisition synchronously. The laser rangefinder sends a trigger pulse signal when the error of the distance to the top of the canopy is less than ±0.2m. S13: The RTK-GPS module on the drone records the longitude and latitude coordinates in real time. The barometric altimeter updates the altitude value every 0.5 seconds. The Beidou timing chip marks the acquisition timestamp with 50ms accuracy. S14: The multispectral camera completes the spectral reflectance acquisition of five bands within 20ms after being triggered. At the same time, the ion mobility spectrometer obtains the concentration values ​​of the three elements C, N, and P within a 150ms detection cycle. The two are synchronized at the hardware level through the FPGA chip. S15: Encapsulate the spectral reflectance data, element concentration values, geographic location data and timestamp according to a unified time reference to form an original multi-source data set.

3. The method for intelligent analysis of spatiotemporal changes in forest plant stoichiometric characteristics according to claim 1, characterized in that: The S2 specifically includes: S21: sort the spectral reflectance and element concentration values ​​in the original multi-source data set formed in S1 in time series, extract the single plant observation sequence at each time point and unify it into a sampling interval of one group every 10 minutes; S22: Using the Lomb-Scargle periodic analysis method, we construct a power spectrum distribution diagram in the frequency domain for data sampled at unequal time intervals, and set 0.0417 Hz as the main periodic reference frequency to identify the periodic signal fluctuation components caused by the day-night temperature difference; S23: The sampling weight factor and temperature covariate are introduced to improve the standard Lomb-Scargle algorithm and obtain the modified weighted power spectrum; S24: Based on the weighted periodic model, the temperature-driven signal drift sequence is reconstructed, and the drift trend curve of each observation point is fitted using the least squares method, and then the trend is removed point by point from the original spectral reflectance and element concentration values ​​to obtain standardized data after temperature drift correction; S25: The drift-corrected spectral reflectance and element concentration values ​​are re-bound to the latitude and longitude coordinates, altitude and timestamp to form a standardized multi-source data set.

4. The method for intelligent analysis of spatiotemporal changes in forest plant stoichiometric characteristics according to claim 3 is characterized in that: The S24 specifically includes: S241: Based on the weighted power spectrum of S23, extract the amplitude coefficients of the cosine term and the sine term at the reference frequency of the corresponding main cycle, and reconstruct the drift signal sequence of each observation point. The reconstruction expression is: ,in, For the Drift estimate for each observation point; is the drift signal amplitude; is the reference angular frequency; For the observation time; is the drift phase angle; S242: Use the least squares method to fit the drift sequence and construct the error square function , whose expression is: ,in, For the observations, representing raw spectral reflectance or element concentration values; S243: Perform signal correction operation on each observation point, remove the estimated drift trend from the original observation value, and obtain the drift-corrected standardized data, which is expressed as: ,in, are normalized values ​​after correction for temperature driven drift.

5. The intelligent analysis method of spatiotemporal variation of forest plant stoichiometric characteristics according to claim 1 is characterized in that: The S3 specifically includes: S31: Based on the standardized multi-source data obtained in S2, according to the altitude data corresponding to each observation data point, an initial data sequence along the elevation direction is constructed, and linear interpolation is performed on the elevation sequence to make the spatial interval of the data sequence on the elevation axis uniform; S32: According to the latitude and longitude coordinates and timestamp of the location of the observation data point, find the corresponding satellite remote sensing solar radiation intensity data , calculate the solar radiation intensity value of each data point, and use the inverse distance weighted interpolation method to spatially resample the solar radiation intensity to generate a two-dimensional spatial continuous distribution map of solar radiation intensity; S33: Use soil moisture sensor to collect soil moisture value in target area , and perform spatial Kriging interpolation calculation based on the latitude and longitude positions of the soil moisture sampling points to generate a spatial continuous distribution map of soil moisture; S34: Based on the data sequence after interpolation in the elevation direction, the solar radiation intensity data obtained by S32 interpolation and the soil moisture data obtained by S33 interpolation are superimposed to form a three-dimensional spatial data structure, and a dynamic weight identification value is assigned to each spatial grid unit. ; S35: The three-dimensional spatial data structure reconstructed in S34 and its dynamic weight identifier are uniformly bound to the corresponding latitude and longitude coordinates, elevation position and timestamp to form a space-time cube with a dynamic weight identifier.

6. The method for intelligent analysis of spatiotemporal changes in forest plant stoichiometric characteristics according to claim 1, characterized in that: The S4 specifically includes: S41: Based on the space-time cube formed in S3, the terrain slope value of each spatial unit is calculated according to the latitude and longitude coordinates and elevation values ​​of the spatial unit, and the initial mapping relationship between the slope and the dynamic weight identifier is determined in combination with the dynamic weight identifier value corresponding to the spatial unit; S42: Calculate the migration rate of the elements in each spatial unit according to the change sequence of the element concentration value in the spatial unit over time, and weight the migration rate using the dynamic weight identifier to establish a mapping relationship between the migration rate and the dynamic weight identifier; S43: Based on the initial mapping relationship between slope and weight and the mapping relationship between migration rate and weight, a nonlinear mapping function between terrain slope and element migration rate is established by using a nonlinear regression method; S44: construct an adaptive sliding window with each space unit in the space-time cube as the center, and the size of the sliding window is adaptively adjusted according to the dynamic weight identifier of the central unit. The larger the weight identifier value, the smaller the window size; S45: within the adaptive sliding window, respectively counting the average concentrations of carbon, nitrogen and phosphorus, and calculating the carbon, nitrogen and phosphorus concentration ratio corresponding to the central unit of the window; S46: Slide the adaptive sliding window gradually along the time axis and the space axis, differentiate the carbon-nitrogen-phosphorus ratios of adjacent moments and adjacent spatial units, calculate the spatiotemporal gradient change of the carbon-nitrogen-phosphorus ratio corresponding to each spatial unit, and generate a spatiotemporal gradient change data sequence of the carbon-nitrogen-phosphorus ratio.

7. The method for intelligent analysis of spatiotemporal changes in forest plant stoichiometric characteristics according to claim 6, characterized in that: The expression of the nonlinear mapping function in S43 is: ,in, is the predicted value of element migration rate; is the terrain slope value; is the saturation value parameter of the element migration rate, indicating the theoretical maximum value of the element migration rate; is the nonlinear change rate parameter of the curve; is the turning point parameter of the slope response, indicating the slope value corresponding to the transition from low speed to high speed of element migration rate; is the natural logarithm base constant.

8. The method for intelligent analysis of spatiotemporal changes in forest plant stoichiometric characteristics according to claim 1, characterized in that: The S5 specifically includes: S51: sorting the spatiotemporal gradient change data sequence of the carbon-nitrogen-phosphorus ratio calculated in S4, extracting the gradient change value of each spatial unit in the continuous time series, and forming a time series feature sequence; S52: The Daubechies wavelet basis function is used to perform three-layer discrete wavelet decomposition on the gradient change time series of each spatial unit, decomposing the original sequence into low-frequency approximate components and high-frequency detail components, where the low-frequency component represents long-term trend information and the high-frequency component represents mutation fluctuation information; S53: retain the low-frequency approximate component, remove the noise interference term in the high-frequency detail component, and restore the long-term trend characteristic curve through inverse wavelet transform to obtain the long-term trend component corresponding to each spatial unit; S54: for the high-frequency detail component, performing correlation analysis between spatially adjacent units, and constructing a distribution map of the high-frequency component in three-dimensional space; S55: Using spatial coordinates as variables, spatial Kriging interpolation is performed on the high-frequency detail components to estimate the propagation path and intensity of sudden fluctuations in space, and to construct a spatially continuous short-term fluctuation component field.

9. The method for intelligent analysis of spatiotemporal changes in forest plant stoichiometric characteristics according to claim 8, characterized in that: The S55 specifically includes: S551: Based on the high-frequency detail component distribution map of the spatial unit obtained in S54, the spatial unit with a significant fluctuation peak is taken as the initial source point to determine the starting position of the spatial propagation of the mutation fluctuation; S552: Taking the source point as the spatial interpolation center, define the interpolation influence radius, extract the high-frequency detail component values ​​of adjacent spatial units within the radius, and use the variation function method to calculate the structural correlation degree of the high-frequency detail component in space, and determine the dominant direction of spatial propagation of the high-frequency detail component; S553: ​​Perform spatial Kriging interpolation step by step along the determined spatial propagation advantage direction to obtain a continuous prediction field of the spatial distribution of high-frequency detail components, and identify the fluctuation propagation path according to the spatial variation trend of the prediction field value; S554: Calculate the wave propagation intensity based on the spatial interpolation result of the high-frequency detail component, determine the numerical difference of the high-frequency detail component between adjacent spatial units on the propagation path, calculate the fluctuation intensity change gradient, and clearly indicate the degree of change of the fluctuation intensity on the propagation path; S555: Integrate the wave propagation path and the wave propagation intensity gradient change to form a complete mutation wave space propagation path map and intensity change field, and output the spatially continuous short-term wave component.

10. The method for intelligent analysis of spatiotemporal changes in forest plant stoichiometric characteristics according to claim 1, characterized in that: The S6 specifically includes: S61: Based on the long-term trend component obtained by decoupling in S5, the initial three-dimensional surface of the long-term trend characteristics is constructed by using the three-dimensional visualization method, with the time axis and altitude axis as the horizontal and vertical coordinates and the long-term trend component value as the vertical coordinate; S62: Using the short-term fluctuation component obtained in S5, with the same three-dimensional coordinate axis system, by superimposing the short-term fluctuation values, a complete three-dimensional map of the trend component and the fluctuation component is formed; S63: Perform spatial derivative analysis on the three-dimensional map after the trend component and the fluctuation component are superimposed, calculate the derivative value of the element flux in the map with the change of altitude, determine the spatial position where the derivative in the altitude direction has a significant mutation, and mark it as the terrain turning point; S64: Taking the determined terrain turning point as the center, analyzing the element flux change trend at the location along the time axis, calculating the flux change rate, and setting a significance threshold, thereby determining the element flux mutation threshold; S65: The determined terrain turning points and the corresponding element flux mutation thresholds are clearly marked on the three-dimensional map and visualized through obvious graphic labels.

Images