Intelligent analysis method for spatiotemporal changes in stoichiometric characteristics of forest plants

Abstract

The present invention relates to the fields of ecological environment monitoring and forest remote sensing technology, and specifically to an intelligent analysis method for the spatiotemporal variation of forest plant stoichiometric characteristics. The method comprises the following steps: S1: obtaining the spectral reflectance and element concentration values of plant leaves in a target area; S2: performing spatiotemporal registration processing; S3: constructing a three-dimensional spatiotemporal matrix to generate a spatiotemporal cube with dynamic weight identifiers; S4: calculating the spatiotemporal gradient variation of the carbon, nitrogen, and phosphorus ratio within each spatial unit; S5: performing feature decoupling on the spatiotemporal gradient variation; and S6: generating a three-dimensional map of stoichiometric characteristics with dual altitude-time coordinate axes. By constructing a multi-factor-driven three-dimensional spatiotemporal analysis model, the present invention achieves decoupling analysis of long-term trends and short-term fluctuations in forest plant stoichiometric characteristics, accurately identifying element flux mutation thresholds at topographic turning points, and effectively enhancing the dynamic monitoring capabilities of ecological processes.

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 spatiotemporal changes in 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 is becoming 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, due to the complexity of the terrain and the coupling of ecological factors, the observation data often suffer from 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 lack the ability to integrate multi-source data and correct for interference from diurnal temperature differences, making it difficult to accurately analyze carbon, nitrogen, and phosphorus ratios across multiple environmental gradients. Furthermore, interpolation methods for multidimensional environmental factors such as terrain slope, soil moisture, and solar radiation are often simplistic, failing to separate long-term trends from short-term fluctuations. In the absence of complete spatiotemporal alignment and feature decoupling, researchers struggle to identify significant changes in element migration rates within forest ecosystems or quantify element flux mutation thresholds at topographic transitions. To overcome these shortcomings and meet the demands for precise monitoring and early warning of forest ecological processes, an intelligent analysis method for spatiotemporal changes in forest plant stoichiometric characteristics is urgently needed. 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 for spatiotemporal changes in forest plant stoichiometric characteristics includes the following steps:

[0006] S1: Using a multispectral drone equipped with an ion mobility spectrometer, the spectral reflectance, element concentration values, and corresponding latitude and longitude coordinates, altitude, and acquisition timestamps of plant leaves in the target area are simultaneously acquired to form an original multi-source dataset.

[0007] S2: Perform spatiotemporal registration processing on the original multi-source dataset, 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;

[0008] 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, and generate a space-time cube with dynamic weight labels;

[0009] 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 the carbon, nitrogen and phosphorus ratio in each spatial unit is calculated through an adaptive sliding window;

[0010] S5: Decouple the temporal and spatial gradient variations to separate the long-term trend component dominated by geographical factors and the short-term fluctuation component caused by environmental mutations;

[0011] S6: Based on the decoupled trend component and fluctuation component, a three-dimensional map of stoichiometric characteristics with altitude-time dual coordinate axes is generated, and the element flux mutation threshold at the terrain turning point is marked.

[0012] Optionally, the S1 specifically includes:

[0013] S11: A multispectral camera array and an ion mobility spectrometer are integrated into the drone's belly. The multispectral camera includes optical sensors with five wavelengths: 450nm, 550nm, 670nm, 780nm, and 850nm. The ion mobility spectrometer's inlet is equipped with a micro air pump to maintain a 0.5mm distance from the surface of the plant leaf.

[0014] S12: When the drone hovers at a height of 3m above the target plant, the laser rangefinder triggers the multispectral camera and the ion mobility spectrometer to start data acquisition synchronously. The laser rangefinder sends a trigger pulse signal when the distance error detected to the top of the canopy is less than ±0.2m.

[0015] S13: The drone's onboard RTK-GPS module records longitude and latitude coordinates in real time. The barometric altimeter updates the altitude value every 0.5 seconds, and the Beidou timing chip timestamps the acquisition with 50ms accuracy.

[0016] S14: The multispectral camera completes spectral reflectance acquisition in five bands within 20ms after being triggered. Simultaneously, 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 via the FPGA chip.

[0017] S15: Encapsulate the spectral reflectance data, element concentration values, geographic location data, and timestamps according to a unified time base to form an original multi-source data set.

[0018] Optionally, the S2 specifically includes:

[0019] S21: Time series sorting of the spectral reflectance and element concentration values in the original multi-source data set generated in S1, extracting the single plant observation sequence at each time point and unifying it into a sampling interval of 10 minutes;

[0020] 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 diurnal temperature difference.

[0021] S23: The standard Lomb-Scargle algorithm is improved by introducing sampling weight factors and temperature covariates to obtain the modified weighted power spectrum;

[0022] 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. This trend is then removed point by point from the original spectral reflectance and element concentration values to obtain standardized data after temperature drift correction.

[0023] 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 dataset.

[0024] Optionally, the S24 specifically includes:

[0025] S241: Based on the weighted power spectrum of S23, extract the amplitude coefficients of the cosine and sine terms at the reference frequency of the corresponding main cycle, and reconstruct the drift signal sequence of each observation point. The reconstruction expression is:

[0026] ,in, For the Drift estimate of observation points; is the drift signal amplitude; is the reference angular frequency; For the observation time; is the drift phase angle;

[0027] 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;

[0028] S243: Perform signal correction operation on each observation point to remove the estimated drift trend from the original observation value, and obtain the drift-corrected normalized data, which is expressed as: ,in, are normalized values after correction for temperature-driven drift.

[0029] Optionally, the S3 specifically includes:

[0030] 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. ;

[0031] S32: According to the latitude and longitude coordinates and timestamp 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;

[0032] 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;

[0033] S34: Based on the data sequence after elevation interpolation, 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. ;

[0034] 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.

[0035] Optionally, the S4 specifically includes:

[0036] 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 corresponding to the spatial unit;

[0037] S42: Calculate the element migration rate in each spatial unit according to the time-varying sequence of the element concentration values in the spatial unit, and weight the migration rate using the dynamic weight identifier to establish a mapping relationship between the migration rate and the dynamic weight identifier;

[0038] 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 using a nonlinear regression method;

[0039] S44: construct an adaptive sliding window with each spatial unit in the space-time cube as the center. 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.

[0040] 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 center unit of the window;

[0041] 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.

[0042] 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 base constant of natural logarithm.

[0043] Optionally, the S5 specifically includes:

[0044] 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;

[0045] S52: Use Daubechies wavelet basis function 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 components represent long-term trend information and the high-frequency components represent sudden fluctuation information;

[0046] 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;

[0047] 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;

[0048] S55: Using spatial coordinates as variables, perform spatial Kriging interpolation on the high-frequency detail components to estimate the spatial propagation path and intensity of sudden fluctuations, and construct a spatially continuous short-term fluctuation component field.

[0049] Optionally, the S55 specifically includes:

[0050] 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 used as the initial source point to determine the starting position of the spatial propagation of the sudden fluctuation;

[0051] 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 spatial structural correlation of the high-frequency detail components and determine the dominant direction of spatial propagation of the high-frequency detail components;

[0052] 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 based on the spatial variation trend of the prediction field values;

[0053] S554: Calculate the wave propagation intensity based on the spatial interpolation results of the high-frequency detail components, determine the numerical difference of the high-frequency detail components 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;

[0054] S555: Integrate the wave propagation path and the wave propagation intensity gradient change to form a complete sudden change wave spatial propagation path map and intensity change field, and output the spatially continuous short-term wave component.

[0055] Optionally, the S6 specifically includes:

[0056] 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 using a 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;

[0057] S62: Using the short-term fluctuation components obtained in S5, and using 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;

[0058] 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 location where the derivative has a significant mutation in the altitude direction and mark it as the topographic turning point.

[0059] S64: Analyze the element flux change trend at the determined terrain turning point along the time axis, calculate the flux change rate, and set a significance threshold to determine the element flux mutation threshold;

[0060] S65: Clearly mark the determined terrain turning points and the corresponding element flux mutation thresholds on the three-dimensional map, and visualize them through obvious graphic labels.

[0061] Beneficial effects of the present invention:

[0062] This method, by fusing multi-sensor observation data with an improved diurnal temperature drift correction algorithm, constructs a three-dimensional space-time cube and uses a composite algorithm combining wavelet transform and spatial kriging interpolation to effectively integrate the influence of multidimensional ecological factors on the stoichiometric characteristics of forest plants. This solution significantly improves the accuracy and resolution of spatiotemporal data under multi-factor gradients such as elevation, solar radiation, and soil moisture, enabling the simultaneous identification of long-term trends and short-term fluctuations, and providing clear spatiotemporal analysis capabilities for the dynamic evolution of carbon-nitrogen-phosphorus ratios.

[0063] 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

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

[0065] Figure 1 Schematic diagram of an intelligent analysis method according to an embodiment of the present invention;

[0066] Figure 2 Schematic diagram of the calculation process of the spatiotemporal gradient variation according to an embodiment of the present invention. DETAILED DESCRIPTION

[0067] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments. It is also noted that, to provide a more detailed description, the following embodiments are best and preferred embodiments, and those skilled in the art may employ alternative methods for implementing certain known technologies. Furthermore, the accompanying drawings are intended only to provide a more detailed description of the embodiments and are not intended to limit the present invention.

[0068] It should be noted that references in the specification to "one embodiment," "an embodiment," "exemplary embodiments," "some embodiments," etc. indicate that the described embodiments may include specific features, structures, or characteristics, but not necessarily every embodiment will include such specific features, structures, or characteristics. Furthermore, when specific features, structures, or characteristics are described in conjunction with an embodiment, it is within the knowledge of persons skilled in the relevant art to implement such features, structures, or characteristics in conjunction with other embodiments (whether or not explicitly described).

[0069] In general, terms can be understood, at least in part, from their 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.

[0070] like Figure 1-Figure 2 As shown in FIG, the intelligent analysis method for spatiotemporal changes in forest plant stoichiometric characteristics includes the following steps:

[0071] S1: Using a multispectral drone equipped with an ion mobility spectrometer, the spectral reflectance, element concentration values, and corresponding latitude and longitude coordinates, altitude, and acquisition timestamps of plant leaves in the target area are simultaneously acquired to form an original multi-source dataset.

[0072] S2: Perform spatiotemporal registration processing on the original multi-source dataset, 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;

[0073] 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, and generate a space-time cube with dynamic weight labels;

[0074] 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 the carbon, nitrogen and phosphorus ratio in each spatial unit is calculated through an adaptive sliding window;

[0075] S5: Decouple the temporal and spatial gradient variations to separate the long-term trend component dominated by geographical factors and the short-term fluctuation component caused by environmental mutations;

[0076] S6: Based on the decoupled trend component and fluctuation component, a three-dimensional map of stoichiometric characteristics with altitude-time dual coordinate axes is generated, and the element flux mutation threshold at the terrain turning point is marked.

[0077] S1 specifically includes:

[0078] S11: A multispectral camera array and an ion mobility spectrometer are integrated into the drone's belly. The multispectral camera includes optical sensors with five wavelengths: 450nm, 550nm, 670nm, 780nm, and 850nm. The ion mobility spectrometer's inlet is equipped with a micro air pump to maintain a 0.5mm distance from the surface of the plant leaf.

[0079] S12: When the drone hovers at a height of 3m above the target plant, the laser rangefinder triggers the multispectral camera and the ion mobility spectrometer to start data acquisition synchronously. The laser rangefinder sends a trigger pulse signal when the distance error detected to the top of the canopy is less than ±0.2m.

[0080] S13: The drone's onboard RTK-GPS module records longitude and latitude coordinates in real time. The barometric altimeter updates the altitude value every 0.5 seconds, and the Beidou timing chip timestamps the acquisition with 50ms accuracy.

[0081] S14: The multispectral camera completes spectral reflectance acquisition in five bands within 20ms after being triggered. Simultaneously, 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 via the FPGA chip.

[0082] S15: Encapsulate the spectral reflectance data, element concentration values, geographic location data and timestamps according to a unified time base 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.

[0083] S2 specifically includes:

[0084] S21: Time series sorting of the spectral reflectance and element concentration values in the original multi-source data set generated in S1, extracting the single plant observation sequence at each time point and unifying it into a sampling interval of 10 minutes;

[0085] S22: Using the Lomb-Scargle periodic analysis method, we construct a frequency-domain power spectrum distribution diagram for data sampled at unequal time intervals. We set 0.0417 Hz (i.e., a 24-hour period) as the primary periodic reference frequency to identify the periodic signal fluctuation components caused by the diurnal temperature difference.

[0086] S23: Introducing sampling weight factors and temperature covariates to improve the standard Lomb-Scargle algorithm and obtain the modified 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 weather station in the same period to improve the sensitivity and accuracy of the period analysis;

[0087] The standard Lomb-Scargle algorithm is improved as follows:

[0088] S231: Set the temperature change rate covariate to , whose expression is: ,in, For the moment External ambient temperature;

[0089] S232: Based on Calculate the sampling weight factor for each observation point , defined as follows: ,in, For the Sampling time, is the sensitivity control coefficient, is the absolute value of the temperature change rate at that moment;

[0090] S233: Introducing the Lomb-Scargle algorithm, the original observations Weighted, modified weighted power spectrum The calculation formula is:

[0091] ,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.

[0092] 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. This trend is then removed point by point from the original spectral reflectance and element concentration values to obtain standardized data after temperature drift correction.

[0093] 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 dataset, which serves as the input basis for subsequent spatiotemporal analysis steps. The above steps improve the algorithm's modeling accuracy of diurnal temperature-driven errors by introducing temperature weights and covariate correction mechanisms into the Lomb-Scargle periodic algorithm, combined with reverse elimination of drift trends, ensuring the comparability and stability of multi-source observation data in the temporal dimension, and laying a solid data foundation for the subsequent construction of spatial interpolation structures and feature decoupling.

[0094] S24 specifically includes:

[0095] S241: Based on the weighted power spectrum of S23, extract the amplitude coefficients of the cosine and sine terms at the reference frequency of the corresponding main cycle, and reconstruct the drift signal sequence of each observation point. The reconstruction expression is:

[0096] ,in, For the Drift estimate of observation points; is the drift signal amplitude; is the reference angular frequency, corresponding to the diurnal cycle; For the observation time; is the drift phase angle;

[0097] 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;

[0098] S243: Perform signal correction operation on each observation point to remove the estimated drift trend from the original observation value, and obtain the drift-corrected normalized 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.

[0099] S3 specifically includes:

[0100] 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. ;

[0101] S32: According to the latitude and longitude coordinates and timestamp 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;

[0102] 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 between data points in meters; : distance attenuation coefficient, the value range is ; The number of adjacent data points used for interpolation calculation is fixed at 8;

[0103] 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;

[0104] The Kriging interpolation formula is: ,in: For location The soil moisture content after interpolation at For the Soil moisture values measured at the sampling points; 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;

[0105] S34: Based on the data sequence after elevation interpolation, 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 indicator 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 altitude 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;

[0106] 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; 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 foundation for subsequent nonlinear relationship modeling and trend decoupling analysis.

[0107] S4 specifically includes:

[0108] 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 corresponding to the spatial unit;

[0109] The formula for calculating the terrain slope value is: ,in: is the terrain slope value, in radians; is the elevation difference between the target spatial unit and the adjacent spatial unit, in meters; is the horizontal distance between the target space unit and the adjacent space unit, in meters;

[0110] The method for determining the initial slope-weight mapping relationship in S41 is as follows:

[0111] First, 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;

[0112] Then, based on the spline function, the slope value is determined Identify values for independent variables and dynamic weights is the initial nonlinear mapping relationship of the dependent variable. The specific calculation adopts the cubic spline interpolation method, in which the number and position of nodes are determined by the uniform distribution of data point density;

[0113] Finally, the least squares criterion is used to optimize the parameters of the spline function so that the sum of squares of the errors between the predicted dynamic weight identification value and the actual observation value is minimized.

[0114] S42: Calculate the element migration rate in each spatial unit according to the time-varying sequence of the element concentration values in the spatial unit, and weight the migration rate using the dynamic weight identifier to establish a mapping relationship between the migration rate and the dynamic weight identifier;

[0115] The calculation formula for element migration rate 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;

[0116] The method for determining the migration rate-weight mapping relationship in the above S42 is as follows:

[0117] First, calculate the element migration rate data obtained by 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;

[0118] 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 squares 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.

[0119] 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.

[0120] 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 using a nonlinear regression method;

[0121] S44: construct an adaptive sliding window with each spatial unit in the space-time cube as the center. 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.

[0122] 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 center unit of the window;

[0123] S46: Slide the adaptive sliding window gradually along the time axis and the space axis, differentiate the carbon, nitrogen, and phosphorus ratios of adjacent moments and adjacent spatial units, calculate the spatiotemporal gradient change of the carbon, nitrogen, and phosphorus ratio corresponding to each spatial unit, and generate a spatiotemporal gradient change data sequence of the carbon, nitrogen, and phosphorus ratio; the above steps establish an accurate nonlinear mapping relationship between terrain slope and element migration rate by utilizing dynamic weight identification, 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, and phosphorus element ratio changes, and lays a reliable data support for the precise decoupling analysis of subsequent trends and fluctuation characteristics.

[0124] 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; is the saturation value parameter of the element migration rate, indicating the theoretical maximum value of the 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; is the base constant of natural logarithm (approximately 2.71828), dimensionless.

[0125] S5 specifically includes:

[0126] 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;

[0127] S52: Use Daubechies wavelet basis function 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 components represent long-term trend information and the high-frequency components represent sudden fluctuation information;

[0128] The calculation formula of three-layer discrete wavelet decomposition is: ,in, is the time series of the gradient change of the original carbon, chlorine and 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 sudden fluctuation information of different scales; is the number of discrete wavelet decomposition layers, set to 3; is the time series index;

[0129] 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;

[0130] 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;

[0131] 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;

[0132] The spatial correlation analysis formula is: ,in, For space unit With adjacent space units The spatial correlation coefficient of the high-frequency detail components between Space units With space unit In the The value of the high-frequency detail component at each time point; Space units 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;

[0133] S55: Using spatial coordinates as variables, spatial Kriging interpolation is performed on the high-frequency detail components to estimate the spatial propagation path and intensity of the mutation fluctuations, and to construct a spatially continuous short-term fluctuation component field. The above steps combine wavelet transform with spatial Kriging interpolation to perform a 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.

[0134] S55 specifically includes:

[0135] 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 used as the initial source point to determine the starting position of the spatial propagation of the sudden fluctuation;

[0136] 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 spatial structural correlation of the high-frequency detail components and determine the dominant direction of spatial propagation of the high-frequency detail components;

[0137] The variogram method includes the following sub-steps:

[0138] Taking the source point as the center, determine multiple spatial sampling point pairs in different directions;

[0139] Calculate the square mean of the difference between the high-frequency detail components of each pair of spatial units in each direction;

[0140] The direction with the smallest square mean is selected as the dominant direction of spatial propagation of high-frequency detail components.

[0141] 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 based on the spatial variation trend of the prediction field values;

[0142] S554: Calculate the wave propagation intensity based on the spatial interpolation results of the high-frequency detail components, determine the numerical difference of the high-frequency detail components 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;

[0143] S555: The fluctuation propagation path and the gradient change of the fluctuation propagation intensity 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.

[0144] S6 specifically includes:

[0145] 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 using a 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;

[0146] 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's stoichiometric characteristics in the altitude and time dimensions;

[0147] 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 location where the derivative has a significant mutation in the altitude direction and mark it as the topographic turning point.

[0148] The calculation formula of spatial derivative is: ,in, is the derivative of the element flux in the altitude direction, in milligrams per kilogram per meter; It is the trend and fluctuation superposition value of element flux in the three-dimensional map; : altitude, in meters; Time, in hours;

[0149] S64: Analyze the element flux change trend at the determined terrain turning point along the time axis, calculate the flux change rate, and set a significance threshold to determine the element flux mutation threshold;

[0150] The calculation formula for the element flux change rate is: ,in: is the element flux change rate, in milligrams per kilogram per hour; Indicates the altitude corresponding to 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;

[0151] The steps for determining the element flux mutation threshold are as follows:

[0152] S641: At the altitude corresponding to the terrain turning point, collect element flux change data for 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;

[0153] 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 to measure whether the element flux change exceeds the normal fluctuation range;

[0154] 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.

[0155] 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, thereby forming 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.

[0156] The present invention encompasses any alternatives, modifications, equivalents, and solutions that fall within the spirit and scope of the present invention. To provide a thorough understanding of the present invention, specific details are described in detail below in connection with the preferred embodiments of the present invention, but those skilled in the art will be able to fully understand the present invention without these detailed descriptions. Furthermore, to avoid unnecessary confusion regarding the essence of the present invention, well-known methods, processes, procedures, components, and circuits have not been described in detail.

[0157] 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 principles of the present invention. These improvements and modifications should also be regarded as within the scope of protection of the present invention.

Claims

1. An intelligent analysis method for spatiotemporal changes in forest plant stoichiometric characteristics, characterized by: The following steps are involved: S1: Using a multispectral drone equipped with an ion mobility spectrometer, the spectral reflectance, element concentration values, and corresponding latitude and longitude coordinates, altitude, and acquisition timestamps of plant leaves in the target area are simultaneously acquired to form an original multi-source dataset. S2: Perform spatiotemporal registration on the original multi-source dataset, 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; introduce the sampling weight factor and temperature covariate to improve the standard Lomb-Scargle algorithm to 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 weather station in the same period to improve the sensitivity and accuracy of the period analysis; 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, 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 the carbon, nitrogen and phosphorus ratio in each spatial unit is calculated through an adaptive sliding window; S5: Decouple the temporal and spatial gradient variations 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 altitude-time dual coordinate axes 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, characterized in that: Said S1 specifically includes: S11: A multispectral camera array and an ion mobility spectrometer are integrated into the drone's belly. The multispectral camera includes optical sensors with five wavelengths: 450nm, 550nm, 670nm, 780nm, and 850nm. The ion mobility spectrometer's inlet is equipped with a micro air pump to maintain a 0.5mm distance from the surface of the plant leaf. S12: When the drone hovers at a height of 3m above the target plant, the laser rangefinder triggers the multispectral camera and the ion mobility spectrometer to start data acquisition synchronously. The laser rangefinder sends a trigger pulse signal when the distance error detected to the top of the canopy is less than ±0.2m. S13: The drone's onboard RTK-GPS module records longitude and latitude coordinates in real time. The barometric altimeter updates the altitude value every 0.5 seconds, and the Beidou timing chip timestamps the acquisition with 50ms accuracy. S14: The multispectral camera completes spectral reflectance acquisition in five bands within 20ms after being triggered. Simultaneously, 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 via the FPGA chip. S15: Encapsulate the spectral reflectance data, element concentration values, geographic location data, and timestamps according to a unified time base to form an original multi-source data set.

3. The intelligent analysis method for spatiotemporal changes in forest plant stoichiometric characteristics according to claim 1 is characterized in that: The S2 specifically includes: S21: Time series sorting of the spectral reflectance and element concentration values in the original multi-source data set generated in S1, extracting the single plant observation sequence at each time point and unifying it into a sampling interval of 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 diurnal temperature difference. S23: The standard Lomb-Scargle algorithm is improved by introducing sampling weight factors and temperature covariates to 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. This trend is then 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 dataset.

4. The method for intelligent analysis of spatiotemporal changes in forest plant stoichiometric characteristics according to claim 3, characterized in that: The S24 specifically includes: S241: Based on the weighted power spectrum of S23, extract the amplitude coefficients of the cosine and sine terms 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 of observation points; 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 to remove the estimated drift trend from the original observation value, and obtain the drift-corrected normalized data, which is expressed as: ,in, are normalized values after correction for temperature-driven drift.

5. The intelligent analysis method for spatiotemporal changes in 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, construct an initial data sequence along the elevation direction according to the elevation data corresponding to each observation data point, and perform linear interpolation 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 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 elevation interpolation, 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 corresponding to the spatial unit; S42: Calculate the element migration rate in each spatial unit according to the time-varying sequence of the element concentration values in the spatial unit, 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 using a nonlinear regression method; S44: construct an adaptive sliding window with each spatial unit in the space-time cube as the center. 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 center 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 base constant of natural logarithm.

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: Use Daubechies wavelet basis function 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 components represent long-term trend information and the high-frequency components represent sudden 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, perform spatial Kriging interpolation on the high-frequency detail components to estimate the spatial propagation path and intensity of sudden fluctuations, and construct a spatially continuous short-term fluctuation component field.

9. The intelligent analysis method for 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 used as the initial source point to determine the starting position of the spatial propagation of the sudden 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 spatial structural correlation of the high-frequency detail components and determine the dominant direction of spatial propagation of the high-frequency detail components; 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 based on the spatial variation trend of the prediction field values; S554: Calculate the wave propagation intensity based on the spatial interpolation results of the high-frequency detail components, determine the numerical difference of the high-frequency detail components 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 sudden change wave spatial propagation path map and intensity change field, and output the spatially continuous short-term wave component.

10. The intelligent analysis method for 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 using a 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, and using 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 location where the derivative has a significant mutation in the altitude direction and mark it as the topographic turning point. S64: Analyze the element flux change trend at the determined terrain turning point along the time axis, calculate the flux change rate, and set a significance threshold to determine the element flux mutation threshold; S65: Clearly mark the determined terrain turning points and the corresponding element flux mutation thresholds on the three-dimensional map, and visualize them through obvious graphic labels.

Images