Grassland ecosystem function evaluation method and system based on multi-source data

Abstract

This invention relates to the field of image processing technology, and more particularly to a method and system for evaluating grassland ecosystem function based on multi-source data. The method includes the following steps: acquiring remote sensing image data and aerial images of the target area, and performing image geometric registration; extracting vegetation spatial distribution indicators of the target area based on the spectral characteristics of ground features obtained from the geometric registration results; retrieving land surface temperature indicators, and combining the land surface temperature indicator retrieving results with acquired meteorological observation data to construct a land surface energy balance curve; identifying functional trends in each season based on the land surface energy balance curve; dividing functional feature domains according to the spatial clustering results of functional trends in each season, and determining the distribution of functional feature domains in each temporal phase; and evaluating the grassland ecosystem function score based on the distribution of functional feature domains in each temporal phase. Compared with traditional evaluation methods based on single-source data, this invention achieves a more accurate, quantitative, and dynamic evaluation of ecological function.

Description

Technical Field

[0001] This invention relates to the field of image processing technology, and in particular to a method and system for evaluating grassland ecosystem function based on multi-source data. Background Technology

[0002] The assessment of grassland ecosystem function mainly relies on ground surveys and monitoring of single ecological indicators, such as biomass estimation, community composition analysis, or soil physicochemical property determination. Although these methods can reflect the ecological state at a local scale, they are limited by spatial coverage, sampling frequency, and human interference, making it difficult to accurately reveal the spatiotemporal dynamic characteristics of grassland ecosystems at regional or even global scales.

[0003] Researchers have begun to explore the use of satellite remote sensing data, meteorological observation data, and ecological monitoring data to comprehensively assess grassland ecological functions. For example, they monitor vegetation cover changes using NDVI (Normalized Difference Vegetation Index), estimate net primary productivity (NPP) using light energy use efficiency models, and extrapolate grassland water use efficiency by incorporating meteorological factors. However, these methods are mostly based on single-source data or simple coupled models, making it difficult to effectively integrate multi-source information at different spatiotemporal resolutions and observation scales. Especially when facing complex terrain and highly heterogeneous grassland environments, a single data source often fails to accurately represent the structure-function coupling relationship of the ecosystem, leading to significant uncertainty in functional assessment. Summary of the Invention

[0004] Therefore, it is necessary for the present invention to provide a method and system for evaluating grassland ecosystem functions based on multi-source data, in order to solve at least one of the above-mentioned technical problems.

[0005] To achieve the above objectives, a grassland ecosystem function evaluation method based on multi-source data includes the following steps:

[0006] Step S1: Acquire remote sensing image data and aerial images of the target area, and perform geometric registration of the remote sensing image data and aerial images;

[0007] Step S2: Based on the spectral characteristics of ground features obtained from the geometric registration results, extract the spatial distribution indicators of vegetation in the target area;

[0008] Step S3: Invert the surface temperature index based on the spatial distribution index of vegetation, and combine the surface temperature index inversion results with the acquired meteorological observation data to construct the surface energy balance curve;

[0009] Step S4: Based on the time series changes of the surface energy balance curve, identify the functional trends of each season; divide the functional characteristic domains according to the spatial clustering results of the functional trends of each season, and determine the distribution of the functional characteristic domains of each time phase.

[0010] Step S5: Evaluate the grassland ecosystem function score based on the distribution of functional characteristic domains in each time phase.

[0011] This application achieves spatiotemporal integration of multi-source information by simultaneously utilizing remote sensing imagery, aerial photography, meteorological observation data, and vegetation spatial distribution indicators, overcoming the limitations of single data sources in terms of spatial coverage, observation frequency, and scale adaptation. Through precise geometric registration of remote sensing and aerial images, and by pairing vegetation distribution information with environmental factors such as surface temperature and latent heat flux, high-resolution, continuous temporal ecological function data covering the entire target area can be obtained, accurately revealing the spatiotemporal dynamic characteristics of grassland ecosystems at the regional scale. Secondly, this application can simultaneously reflect the physiological activities and cover structure evolution of grassland vegetation. By converting the dynamic trend of latent heat flux into vegetation transpiration activity and combining it with the dynamic trend of vegetation cover, a functional trend matrix is ​​formed. Furthermore, hierarchical clustering and spatial proximity determination are used to divide functional characteristic domains, thereby establishing a structure-function coupled spatial pattern. This approach not only ensures the spatial continuity and rationality of the evaluation results but also enhances the sensitivity and explanatory power for changes in ecosystem function under complex and heterogeneous grassland environments. Furthermore, this application, through comprehensive analysis of functional characteristic domains in both stable and evolving phases, incorporates multi-dimensional indicators such as area proportion, spatial continuity, center offset, and area change rate into the functional scoring system, achieving quantitative and dynamic ecological function evaluation. Compared to traditional evaluation methods based on single-source data, this application has significant advantages in spatiotemporal accuracy, dynamic response capability, and regional applicability, contributing to supporting grassland ecological management, ecological restoration assessment, and regional ecological function classification decision-making.

[0012] Optionally, this specification also provides a grassland ecosystem function evaluation system based on multi-source data, used to perform the grassland ecosystem function evaluation method based on multi-source data as described above. The grassland ecosystem function evaluation system based on multi-source data includes:

[0013] The image registration module is used to acquire remote sensing image data and aerial images of the target area, and to perform geometric registration of the remote sensing image data and aerial images.

[0014] The spatial distribution extraction module is used to extract vegetation spatial distribution indicators of the target area based on the spectral characteristics of ground features from the geometric registration results.

[0015] The energy balance curve construction module is used to invert the surface temperature index based on the vegetation spatial distribution index, and combine the surface temperature index inversion results with the acquired meteorological observation data to construct the surface energy balance curve.

[0016] The functional feature domain segmentation module is used to identify the functional trends of each season based on the time series changes of the surface energy balance curve; to segment the functional feature domains according to the spatial clustering results of the functional trends of each season; and to determine the distribution of the functional feature domains in each time phase.

[0017] The functional score assessment module is used to evaluate the grassland ecosystem functional score based on the distribution of functional characteristic domains in each time phase.

[0018] This invention discloses a grassland ecosystem function evaluation system based on multi-source data. This system can implement any of the grassland ecosystem function evaluation methods based on multi-source data of this invention. It serves as a medium for coordinating the operations and signal transmission between various modules to complete the grassland ecosystem function evaluation method based on multi-source data. The modules within the system cooperate with each other, thereby achieving quantitative and dynamic ecological function evaluation. Compared with traditional evaluation methods based on single-source data, this application helps support grassland ecological management, ecological restoration assessment, and regional ecological function classification decision-making. Attached Figure Description

[0019] Other features, objects, and advantages of the invention will become more apparent from the following detailed description of non-limiting embodiments with reference to the accompanying drawings:

[0020] Figure 1 This is a schematic diagram of the steps in the grassland ecosystem function evaluation method based on multi-source data of the present invention;

[0021] Figure 2 This is a schematic diagram of a vegetation distribution unit in an embodiment of the present invention;

[0022] Figure 3 This is a schematic diagram of the surface energy balance curve in an embodiment of the present invention;

[0023] The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0024] The technical method of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0025] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0026] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0027] To achieve the above objectives, please refer to Figures 1 to 3 This invention provides a method for evaluating grassland ecosystem function based on multi-source data, the method comprising the following steps:

[0028] Step S1: Acquire remote sensing image data and aerial images of the target area, and perform geometric registration of the remote sensing image data and aerial images;

[0029] In this embodiment, satellite remote sensing imagery with a spatial resolution of 0.8m is acquired through a timed scheduling method, while low-altitude aerial photography is performed by a UAV equipped with visible light and near-infrared imaging units. After acquisition, an image geometric registration method based on feature point matching is used to establish a control point set by detecting road intersections, riverbank polylines, and right-angled boundary shapes of facilities in the imagery. Subsequently, a geometric transformation matrix is ​​formed by calculating the pixel offset of the control points. When there are fewer than 5 ground calibration points, the degree of gradient of ground feature boundaries is also checked, and adjustments are made according to the clarity of vegetation edges. Next, the remote sensing imagery and the aerial photographic imagery are aligned according to timestamps, and a time tolerance window of approximately 0.5 seconds is defined to ensure image synchronization. The final output is aligned image data containing a unified coordinate frame and spatial reference, which is used for subsequent vegetation information extraction to form a joint image data sequence with traceable time labels.

[0030] Step S2: Based on the spectral characteristics of ground features obtained from the geometric registration results, extract the spatial distribution indicators of vegetation in the target area;

[0031] In a further embodiment, by detecting the gradient of vegetation reflectance intensity across different wavelengths, the response differences in the near-infrared and red light bands are calculated. During implementation, all pixels in the image are spatially averaged using a 5×5 window to suppress local noise, and edge attenuation parameters are recorded to determine vegetation continuity. Subsequently, the spectral response curves are compared with a previously collected standard plot spectral library, which consists of 50 stable vegetation plots and includes slope parameters in the near-infrared band, red light absorption depth, and seasonal reflectance shift. By calculating the cumulative difference between each pixel and the standard curve, the vegetation coverage area is identified, and vegetation spatial units are generated to obtain a vegetation spatial distribution index matrix. This matrix contains information such as pixel coverage percentage, patch area ratio, and boundary extension direction, which is used for subsequent functional inversion analysis.

[0032] Figure 2 This is a schematic diagram of a vegetation distribution unit in an embodiment of the present invention; as shown below. Figure 2 As shown in the figure, the different colored blocks represent vegetation spatial units. The color differences stem from differences in the spectral response of vegetation species, growth stages, or cover density. Based on the spectral analysis logic of the above embodiments, each color can be distinguished as follows: Green areas (light green, dark green, yellowish-green) represent dominant vegetation cover areas. They have strong near-infrared reflection and significant red light absorption. Their spectral curves closely match those of the "standard vegetation sample area," representing vegetation patches with good continuity and high cover density, such as dense woodlands and well-grown grasslands. Blue-green areas represent wetland vegetation or special vegetation types. Their spectral responses exhibit unique differences in the near-infrared and red light bands (e.g., high water content leads to a slight decrease in near-infrared reflection). Spatially, they are mostly relatively regular blocks, representing vegetation units with specific ecological attributes. The red / orange color scheme belongs to the vegetation cover transition zone or special vegetation type. Its spectral curves differ from those of the standard plots (e.g., lower vegetation density, containing a small amount of non-vegetation components). Spatially, it is mostly marginal patches with weak vegetation continuity. It is an area of ​​sparse forest, shrubland, or vegetation restoration.

[0033] Step S3: Invert the surface temperature index based on the spatial distribution index of vegetation, and combine the surface temperature index inversion results with the acquired meteorological observation data to construct the surface energy balance curve;

[0034] In a further embodiment, temperature estimation weights are formed based on the thermal radiation intensity of each pixel in the aerial image and the vegetation cover level, and then the surface temperature index is calculated. This is then integrated with relative humidity, wind speed, and solar radiation flux collected from meteorological observation stations. Meteorological data are recorded at 10-minute intervals, with a 2-minute time deviation allowable threshold set; if exceeded, linear time interpolation is performed. After data integration, a surface energy balance curve is constructed. This curve can include components such as net radiation, sensible heat flux, and latent heat flux. By detecting the fluctuation amplitude and peak shift of these components over time, vegetation transpiration and energy transfer efficiency can be characterized. Finally, a set of time-stamped energy balance curves is output, providing a basis for seasonal functional characteristic analysis.

[0035] Step S4: Based on the time series changes of the surface energy balance curve, identify the functional trends of each season; divide the functional characteristic domains according to the spatial clustering results of the functional trends of each season, and determine the distribution of the functional characteristic domains of each time phase.

[0036] In a further embodiment, the slope change of the surface energy balance curve during the morning warming phase can be detected, the fluctuation range of the daily average sensible heat flux can be calculated, and the peak time of the latent heat flux can be determined. If the time deviation exceeds 30 minutes for three consecutive days, the functional performance of the area is considered to have entered a non-steady-state period. Subsequently, the annual energy curve is segmented by season, and any two characteristic parameters from flux peak, energy dissipation width, vegetation transpiration activity, vegetation cover, and thermal response delay are extracted to form a seasonal functional trend feature set. Spatial clustering is performed on this feature set, and functional feature aggregation areas are formed by judging the cumulative difference in parameter values ​​of neighboring pixels, thereby dividing the grassland ecological function feature domain. Finally, the location, boundary, and characteristics of each functional domain under different seasons are labeled with temporal sequence numbers to form a spatially oriented functional domain distribution result.

[0037] Step S5: Evaluate the grassland ecosystem function score based on the distribution of functional characteristic domains in each time phase.

[0038] In a further embodiment, the area proportion and spatial continuity of functional characteristic domains can be used as the basic score, and the area change rate and the center offset rate obtained based on the center position offset can be used as deduction factors to evaluate the grassland ecosystem function score of the target area. This score is quantified in the range of 0 to 100. If the grassland ecosystem function score of any functional characteristic domain is less than 30% of the average grassland ecosystem function score of the target area, an early warning measure is initiated through the grassland ecosystem monitoring platform.

[0039] In another embodiment, the seasonal mean of latent heat flux, vegetation cover ratio, and energy transfer stability within each functional domain can be detected and compared with historical baselines. If the difference exceeds a preset reference range of 0.12, the ecological performance is considered to have deviated. Simultaneously, the shift in peak surface temperature is detected to assess the impact of changes in hydrothermal conditions on grassland growth. Subsequently, weights are assigned according to the temporal persistence, spatial morphological integrity, and energy response continuity of different functional domains, generating an ecosystem function score quantified from 0 to 100. Upon output, the score is bound to the corresponding time phase identifier, forming a grassland ecosystem function evaluation result that can be used for cross-seasonal comparisons.

[0040] Optionally, performing image geometric registration in step S1 includes:

[0041] Detect the edge feature differences between any image frame in remote sensing image data and aerial images, and select control pixels based on the edge feature differences;

[0042] In this embodiment, edge feature detection is performed on any image frame and aerial image from the remote sensing image data. During implementation, the Canny operator is applied to extract edge contours from the image frame and aerial image respectively, with parameters set to a low threshold of 0.1 and a high threshold of 0.3. A 3×3 Gaussian filter is then used for noise smoothing to eliminate local brightness interference. After detection, the edge pixels are combined with texture gradients. By calculating the amplitude and direction consistency of grayscale changes within a 3×3 neighborhood window, the edge feature response value of each pixel is obtained. When the response value exceeds the threshold of 0.25, it is included in the control pixel candidate set, providing a basic data source for subsequent similarity matching.

[0043] Calculate the similarity matrix of control pixels in the image frame and the aerial image respectively to obtain the correspondence of control pixels;

[0044] In this embodiment, a feature vector is constructed based on the candidate control pixels. This vector includes three parameters: pixel grayscale value, texture orientation gradient, and neighborhood uniformity, used to accurately characterize the local structure of the pixel. Candidate control points in the image frame and aerial image are paired and a similarity matrix is ​​calculated. The similarity is obtained by weighted summation of grayscale difference, orientation difference, and neighborhood texture deviation, with weights set to 0.5 for grayscale difference, 0.3 for orientation difference, and 0.2 for texture deviation. If the similarity of a control pixel pair is greater than 0.87, it is considered a valid correspondence, and the coordinate information of the correspondence is recorded in the correspondence matrix for subsequent calculation of spatial mapping parameters.

[0045] Based on the coordinate deviation in the correspondence of control pixels, the spatial mapping parameters between the image frame and the aerial image are determined respectively, and geometric resampling is performed on the pixel gray values ​​of the image frame and the aerial image according to the spatial mapping parameters.

[0046] In a further embodiment, based on the control pixel correspondence matrix, spatial mapping parameters for the image frame and the aerial image are calculated, including translation in the longitudinal and lateral directions, rotation angle, and linear scaling factor. The initial value of the scaling factor is set to 1, and the rotation angle is iteratively adjusted in steps of 0.05° until the reprojection error of the control points is minimized, with an error tolerance set to 0.2 pixels. After obtaining the spatial mapping parameters, geometric resampling processing is performed on the grayscale values ​​of corresponding pixels in the image frame and the aerial image. A bilinear interpolation method can be used, with the interpolation neighborhood set to 2×2 pixels to maintain edge continuity and texture integrity. After resampling, a spatially unified geometrically registered image frame is generated, providing basic data for subsequent extraction of vegetation spatial distribution indicators.

[0047] The geometric resampling results are integrated in time according to the timestamp, and the time integration result is used as the geometric registration result.

[0048] In a further embodiment, the geometric resampling results are sequentially integrated according to the timestamps of the image frames and the aerial images. During the integration process, a time deviation tolerance of ±1 second is set, and missing frames are filled using linear interpolation to ensure the integrity of the time series. The integrated results form a geometrically registered image sequence, in which each image frame and the aerial image correspond to a unified spatial coordinate reference system, including pixel grayscale values ​​and edge feature response information.

[0049] Optionally, selecting control pixels includes:

[0050] Edge detection and texture gradient calculation are performed on any image frame and aerial image in the remote sensing image data, respectively. The edge detection results and texture gradient information of the image frame and aerial image are then normalized and fused according to pixel position to obtain feature response pairs.

[0051] In this embodiment, edge detection based on the Canny operator is performed on remote sensing image frames and aerial images respectively, with the detection threshold range set to 0.05 to 0.18, in order to extract clear ground feature boundary lines. At the same time, the Sobel operator is used to calculate the texture gradient direction and gradient magnitude, and the edge detection results and texture gradient results are linearly normalized and fused according to pixel position to generate feature response pairs. This makes the fused response contain both ground feature outline information and texture detail sensitivity, laying the foundation for subsequent control point selection.

[0052] Based on the texture repetition of each pixel in the feature response pair, regions in image frames and aerial images with texture repetition exceeding a preset texture repetition threshold are removed.

[0053] In a further embodiment, to avoid feature ambiguity caused by repetitive surface textures, texture redundancy analysis is performed on the neighborhood window of each pixel. The texture redundancy index R is calculated by statistically analyzing the local gradient direction variance within a 3×3 pixel neighborhood. If R is greater than a preset texture redundancy threshold of 0.62, the region is considered to have a periodic structure (such as grassland stripes, field ridges, or building roof tiles), and thus the region is removed from the matching range. This process effectively reduces interference from redundant texture regions, enhances the spatial representativeness of feature distribution, and provides a stable basis for the uniform selection of control points.

[0054] The region removal results are divided into grids, and the cell with the highest response value in each grid is selected as a candidate control point.

[0055] In a further embodiment, after region culling, to ensure a balanced distribution of control pixels across the entire space, the remaining feature regions in the image frame and aerial image are divided into regular grid units, each covering 32×32 pixels. Within each grid, the pixel with the highest response value is selected as a candidate control point based on the intensity of the fused feature response value. This selection method avoids the phenomenon of dense control points in local areas due to high-frequency texture aggregation, ensuring the spatial uniformity of control points in the distribution of ground features.

[0056] Calculate the local grayscale values ​​of candidate control points at different resolutions in remote sensing image data and aerial images, and select candidate control points whose local grayscale value standard deviation is less than the average grayscale value of the image as control pixels.

[0057] In a further embodiment, different resolutions (0.5m, 1m, and 2m) are extracted from the remote sensing image data and aerial images corresponding to the candidate control points, and a 5×5 local grayscale window is formed centered on the candidate point. The standard deviation of the grayscale value of each window is then calculated. With the corresponding image average gray level .when Less than When the ratio threshold is set to ×0.12, the candidate point is determined to have stable brightness under multi-scale conditions and be unaffected by resolution changes; if this condition is not met, the candidate point is removed.

[0058] Optionally, step S3, which involves retrieving surface temperature indicators based on vegetation spatial distribution indicators, includes:

[0059] By combining the infrared band information of remote sensing image data in the geometric registration results with vegetation spatial distribution indicators, the surface thermal radiation intensity of each vegetation distribution unit in the vegetation spatial distribution indicators is estimated.

[0060] In this embodiment, the infrared band information of the remote sensing image data is called in the completed geometric registration result and superimposed with the vegetation spatial distribution index generated based on NDVI. In this implementation, NDVI is constructed from the difference between the near-infrared reflectance and red-light reflectance after normalization processing. After pixel-level calculation, it is divided into several vegetation distribution units, and each unit covers approximately 10×10 pixels. When the incident radiation luminance value of the infrared band is within the range of 0.35 to 0.58, the system accumulates the thermal radiation intensity unit by unit according to the spatial position after pixel geometric correction and the vegetation index level, forming a primary radiation intensity matrix. Subsequently, at the same timestamp, a time-synchronized average luminance correction is performed on the matrix, enabling the vegetation distribution unit to present the corresponding surface thermal radiation intensity by virtue of its infrared response ability. This method can avoid the distortion of vegetation energy response caused by shadow coverage or local reflection anomalies.

[0061] Calculate the spectral radiance parameters of each vegetation distribution unit based on the surface thermal radiation intensity, and use the spectral radiance parameters to correct the surface thermal radiation intensity;

[0062] In a further embodiment, the spectral radiance parameters of each vegetation distribution unit are further calculated based on the surface thermal radiation intensity. This embodiment adopts a multi-band comparison method, calculating the difference in reflection intensity between the main infrared band (central wavelength 10.8 ) and the near-infrared auxiliary band (central wavelength 8.7 ). At the same time, referring to the typical spectral response characteristics of leaf epidermal roughness and water content in the on-site historical sample library, the roughness index and the water content index are introduced into the calculation, where they vary within the ranges of 0.2 < r < 0.6 and 0.1 < w < 0.4. By comparing the feedback deviations of the two bands and fitting the reflection-radiation curves of isomorphic plants in the sample library, the spectral radiance parameters of the vegetation unit are obtained. This parameter is written back to the aforementioned surface thermal radiation intensity, performing a proportional correction on the primary radiation value, so as to adjust the thermal response deviation caused by differences in leaf direction and structural density, forming the secondary surface thermal radiation intensity. This process ensures that the micro-structure and moisture changes on the vegetation surface have the ability to be quantitatively corrected.

[0063] Convert the corrected surface thermal radiation intensity into surface temperature, and perform a spatial weighted average on the surface temperature at the same timestamp to obtain the surface temperature index. <00>

[0064] In a further embodiment, the corrected surface thermal radiation intensity is mapped to temperature using a radiation conversion formula. An atmospheric interference compensation coefficient k (ranging from 0.03 to 0.07) is introduced during the mapping process to offset the slight energy shift caused by canopy reflectivity. After the conversion, a surface temperature raster layer is obtained, and a spatially weighted average is performed at the same timestamp. Spatial weighting uses the vegetation density index as a weighting factor; higher density results in a greater weight, reflecting the thermal inertia effect of vegetation. By performing a sliding weighted average within a 9×9 pixel area of ​​a local window, local noise can be balanced while preserving the macroscopic thermal trend, thus obtaining the surface temperature index.

[0065] Optionally, constructing the surface energy balance curve in step S3 includes:

[0066] Acquire meteorological observation data; synchronize the surface temperature index with the meteorological observation data at the same timestamp, and pair them according to the spatial location corresponding to the vegetation distribution unit to obtain a multi-source data set of the unit;

[0067] In this embodiment, the meteorological observation data comes from ground-based automatic weather stations, and the observed elements include wind speed (2.1–3.8 m / s), air humidity (60%–75%), and radiation flux density (450–610). The data included surface temperature and near-surface air temperature (18.5–22.3℃). Meteorological stations were deployed 200m upwind of the vegetated area, using a 1-minute sampling period and recording the data in a timestamp format in the data buffer. Subsequently, using the resampling position of the vegetation distribution unit in the geometric registration coordinate system as a reference, surface temperature indicators and meteorological observation data were matched at the same timestamp. A spatial neighborhood threshold (radius 15m) was introduced during the matching process to automatically retrieve the corresponding meteorological observation point. If multiple sampling records existed within the radius, the strategy of minimizing the absolute value of the timestamp difference was used to lock in the optimal matching record. Finally, a multi-source data set containing surface temperature, wind speed, air humidity, radiation flux density, and air temperature was formed, providing stable input for subsequent thermodynamic quantity calculations.

[0068] Input the multi-source data set of the unit into the preset energy budget equation to calculate the latent heat flux of each vegetation distribution unit;

[0069] In a further embodiment, the energy budget equation is derived from a previously calibrated semi-empirical energy balance structure to express the energy relationship between solar radiation absorption, sensible heat exchange, latent heat evaporation, and surface heat storage. Meteorological humidity is introduced as a constraint on the water vapor gradient, and wind speed is used as a correction factor for near-surface turbulence intensity; specifically, when the wind speed is greater than 2.4 m / s, the turbulence correction factor is 0.12, otherwise it is 0.08. To avoid the model appearing out of thin air, this structure is trained based on an international agricultural meteorological observation sample database and includes an input layer (near-surface air temperature...). ), air humidity (%), wind speed (m / s), surface temperature ( The structure consists of five elements: a radiative flux density, an intermediate coupling layer (coupling air temperature and humidity gradient, surface temperature difference, and radiative response), and an output layer (latent heat flux), with weights derived from historical sample fitting. By substituting each element of the multi-source dataset into this structural formula, the latent heat flux of the vegetation unit can be obtained, ranging from 50 to 120. The latent heat fluxes described above accurately reflect changes in vegetation transpiration intensity and leaf area heat dissipation capacity.

[0070] It is worth noting that the semi-empirical energy balance structure is constructed based on the general principle of surface energy conservation and calibrated using historical meteorological observation samples from the region. Its theoretical basis comes from the energy conservation equation, namely: the net radiative energy absorbed by the surface is distributed among sensible heat flux, latent heat flux, and surface heat storage. Due to evapotranspiration in vegetation-covered areas, this structure adds an evapotranspiration coupling term to the general framework. Historical sample data are selected from continuous observation sequences of meteorological environments with similar vegetation. Parameter regression is used to empirically fit factors such as the influence of evapotranspiration, wind-driven convection, and the influence of humidity gradient on latent heat direction, thereby obtaining semi-empirical parameters adapted to the seasonal characteristics of the region. This semi-empirical energy balance structure can be expressed as: ;in This represents the net radiation flux at the Earth's surface. This is the sensible heat flux (where this term is a semi-empirical convection term and requires calibration to determine the empirical coefficient). For latent heat flux, This refers to surface heat storage.

[0071] The temporal variation of latent heat flux in each vegetation distribution unit is plotted as a surface energy balance curve.

[0072] In this embodiment, the sampling interval (300s) is used as the horizontal axis increment during the plotting process, and latent heat flux is used as the original vertical axis reference. Radiation flux density is introduced as an auxiliary background quantity for bias verification. Before plotting, the latent heat flux sequence is subjected to sliding smoothing, with the smoothing window width set to 3 sampling periods to reduce thermal shifts caused by local meteorological disturbances. Subsequently, the latent heat values ​​at different time points are plotted as a broken line, and key inflection points are marked according to the radiation background to show the changes in the equilibrium state of vegetation energy absorption-evaporation release process. By observing the curves, typical segments such as enhanced daytime transpiration, stabilization in the evening, and slow release of residual heat at night can be identified, enabling the energy budget process of vegetation ecological units to be presented in a visual way, providing a basic interpretation basis for subsequent ecological performance assessment and drought stress detection.

[0073] Optionally, identifying the functional trends for each season in step S4 includes:

[0074] The spatial distribution indicators of vegetation and the surface energy balance curve at the same time stamp are spatiotemporally synchronized according to the vegetation distribution unit, and the spatiotemporal synchronization results are used as the data source to be analyzed.

[0075] In this embodiment, the remote sensing observation period is divided into several timestamps. When constructing the basic data for analysis, the vegetation spatial distribution index is first obtained based on the registered remote sensing images and aerial data from previous steps. At each timestamp, the vegetation distribution index and the surface energy balance curve are aligned according to the spatial grid boundary of the vegetation distribution unit (e.g., 50m × 50m unit resolution). To ensure the consistency of data synchronization, a maximum time offset threshold is set for the observation records of each timestamp. For intervals ≤10 minutes, linear interpolation is used to correct for any offsets, ensuring that vegetation spatial distribution indicators and latent heat flux curves are spatiotemporally integrated on the same vegetation distribution unit. The integrated result is then used as the data source for analysis. Thus, the data source for analysis forms a three-element data structure of "vegetation cover characteristic value + latent heat flux curve value + time series index," without introducing external data, forming a closed and traceable data foundation.

[0076] Based on the data sources to be analyzed, determine the dynamic trends of undercurrent heat flux and vegetation cover in each season;

[0077] In a further embodiment, based on the aforementioned data source to be analyzed, the annual timeline is divided into multiple seasonal windows according to the general meteorological seasonal division rules (e.g., three months per season). The time series changes in latent heat flux are extracted within each seasonal window to depict the seasonal dynamic trend. This trend is calculated using a sliding window (7-day window width, 3-day step size) to measure the increasing or decreasing tendency of the flux. Simultaneously, the spatial mean change rate of vegetation cover within the same window is calculated to reflect the direction of seasonal change in vegetation density. Subsequently, to improve the comparability of trends, zero-mean standardization is applied to both the latent heat flux trend value and the vegetation cover change trend, ensuring that each trend falls within the [-1,1] interval to avoid dimensional differences. This yields dynamic trend arrays of latent heat flux and vegetation cover for each season, providing a basis for subsequent transpiration intensity conversion.

[0078] The dynamic trend of latent heat flux in each season is converted into the vegetation transpiration activity trend by the unit area of ​​each vegetation distribution unit. The vegetation transpiration activity trend and the vegetation cover dynamic trend in each season are then standardized and aggregated as the functional trend of each season.

[0079] In a further embodiment, after obtaining the dynamic trend of latent heat flux over time, in order to reflect the transpiration contribution of vegetation per unit area to surface energy, the unit area coefficient of the corresponding vegetation distribution unit (e.g., 2500) is used. In the conversion process of latent heat dynamic trends (corresponding to 50m×50m), transpiration activity trend values ​​are formed. This conversion uses an energy density normalization operation, that is, dividing the latent heat dynamic trend value by the unit area and multiplying it by a seasonal scale correction factor (set between 0.85 and 0.92 to reduce the interference of abnormal weather on the trend), thus obtaining a trend curve that more closely reflects actual transpiration feedback. At the same time, the vegetation cover dynamic trend is kept consistent in numerical range through the same standardization operation to avoid comparison bias. Finally, the vegetation transpiration activity trend converted per unit area and the standardized vegetation cover dynamic trend under the same season are aggregated in the form of two-dimensional trend pairs to form the functional trend data structure for each season, which is used for subsequent functional feature domain clustering and ecological assessment.

[0080] Optionally, the functional feature domain division in step S4 includes:

[0081] Based on the functional trends of vegetation transpiration activity in each season, the first feature dimension is used as the first feature dimension, and the dynamic trend of vegetation coverage is used as the second feature dimension.

[0082] Calculate the trend similarity of the first feature dimension and the second feature dimension between each vegetation distribution unit, and perform hierarchical clustering based on the trend similarity, merging distribution unit clusters with similarity higher than the preset similarity threshold layer by layer;

[0083] In this embodiment, a two-dimensional trend feature coordinate system is constructed. This system uses the vegetation transpiration activity trend as the vertical feature axis and the vegetation cover change trend as the horizontal feature axis. The trend feature matrix of the distribution units is then established using this two-dimensional result. The matrix uses rows to represent spatial unit numbers, and columns to correspond to the dynamic trend values ​​of transpiration activity and vegetation cover change, respectively. To avoid dimensional bias caused by differences in trend scales, each trend value is standardized before matrix construction to ensure it falls within the range [-1, 1]. For the standardized matrix, the difference between each unit on the two feature axes is obtained by calculating the trend difference between each unit. The difference is expressed using discretized weight parameters. =0.6 imparts a higher proportion of transpiration activity, to =0.4 assigns a low weight to the coverage trend. This weighting is based on empirical analysis that latent heat flux plays a more dominant role in vegetation transpiration in the preceding energy balance structure, ensuring that the trend evaluation constitutes a logical closure. After the difference calculation is completed, the difference values ​​are mapped to similarity scores; the higher the similarity, the more consistent the trend performance.

[0084] Furthermore, the hierarchical aggregation structure used in this embodiment is derived from a commonly used agglomerative hierarchical aggregation strategy. In its initial stage, this structure treats each vegetation distribution unit as an independent cluster. Then, using a trend similarity matrix as input data, it recursively merges the most similar clusters. This structure uses a dendrogram-like association method to record the merging process, and its core function is to dynamically evaluate the trend differences between different clusters and construct an aggregation relationship tree from the bottom up. Simultaneously, to control the merging sensitivity, a preset trend similarity threshold is introduced. The similarity is set to 0.78. If the trend similarity between clusters to be merged is below this value, further merging is not allowed to avoid unreasonable integration. To ensure spatial continuity and avoid clustering of distant units due to aggregation based solely on eigenvalues, the distance matrix is ​​called and the spatial proximity score is calculated for each aggregation operation. ,like If the value is less than the constraint threshold of 0.65, the current round of merging is cancelled, thereby ensuring that the spatial regions corresponding to the merged clusters are distributed in a continuous pattern.

[0085] Based on the merging results of the distribution unit clusters in the first feature dimension and the second feature dimension, vegetation distribution units that simultaneously satisfy the conditions of spatial proximity and show a unidirectional change relationship in both the first and second feature dimensions are divided into functional feature domains.

[0086] In a further embodiment, when the tree-like aggregation structure is completed and stops at the similarity threshold... When the system reaches a certain point, it obtains the final distribution unit cluster division. Subsequently, this division result is compared and verified with the trend aggregation result of the second feature dimension to determine whether there is a co-directional change in the two aggregation results. Co-directional change here refers to the simultaneous positive increase or simultaneous decline of transpiration activity and cover within the same period, determined by the sign parameter. =+1 or =-1 indicates directionality to ensure consistency in trend interpretation. If a group of units is identified in the control group that simultaneously satisfies spatial proximity, similar trends, and unidirectional change, this group of units is classified into the same functional domain to form a spatial region that reflects the coupling between vegetation physiological activity and surface cover evolution. During process execution, density maintenance parameters are also used... This is used to limit the minimum area contained in the feature domain, ensuring that excessively small clustering results do not generate noise. In this implementation example... The value is set to 0.12, meaning that if the area of ​​the formed candidate aggregation region is smaller than this ratio, it will automatically revert to the previous aggregation state to avoid the generation of fragmented functional domains. The functional feature domains obtained through the above process can be used for subsequent surface water and heat exchange analysis, vegetation stress response studies, or regional ecological function classification assessments, forming a closed and complete data processing logic chain.

[0087] Of particular importance are the methods for determining spatial proximity, which include:

[0088] Based on the geographic coordinates of each vegetation distribution unit, the Euclidean distance of all vegetation distribution units in the cluster merging result is calculated; if the calculated Euclidean distance of any vegetation distribution unit pair is less than the preset spatial proximity threshold, the spatial relationship of the vegetation distribution unit pair is determined to be spatially adjacent.

[0089] In this embodiment, the geographic coordinates of all vegetation distribution units are extracted from the vegetation spatial distribution index. These coordinates can be latitude and longitude or planar coordinates in a projected coordinate system. Then, for each pair of vegetation distribution units within a cluster, their Euclidean distance is calculated to form an Euclidean distance matrix. To improve computational efficiency, the Euclidean distance calculation process can be combined with a spatial indexing structure (such as a quadtree or kd-tree) to achieve fast neighborhood search, which is particularly suitable for areas containing thousands to tens of thousands of vegetation units. Next, the Euclidean distance matrix is ​​compared with a preset spatial proximity threshold. The threshold can be set based on the remote sensing image resolution and the target area scale. For example, if the remote sensing image resolution is 10 meters, the threshold can be... The distance is set to 30 meters to ensure significant spatial physical connections between units. If the Euclidean distance between any pair of vegetation units is less than a threshold... If the unit satisfies the condition of spatial proximity, it will be included in the same cluster in the subsequent functional feature domain division to ensure that the formed functional feature domains are spatially continuous and have clear physical connections, while avoiding the generation of isolated or fragmented regions.

[0090] Optionally, determining the distribution of functional characteristic domains for each time phase in step S4 includes:

[0091] The preset remote sensing observation period is used as the time phase, and the geographic boundary contour of each functional feature domain is extracted in the current time phase based on the spatial position relationship of vegetation distribution units in the functional feature domain.

[0092] In this embodiment, a preset remote sensing observation period is used. As a temporal unit, the latitude and longitude coordinates and number information of each vegetation distribution unit within the functional feature domain are time-stamped and integrated to construct a multi-temporal spatial index table. This index table is structured in tabular form, with rows representing the number of each vegetation unit and columns representing the latitude and longitude coordinates and unit area under different observation time phases, facilitating subsequent geographic boundary contour extraction. Subsequently, the Convex Hull algorithm or... shape( The -shape method generates a contour of the set of unit coordinates within the functional feature domain, resulting in the geographic boundary contour of the functional feature domain for each time phase. The contour is stored in the form of vector polygons, and each polygon node contains x and y coordinates and the corresponding unit number to ensure a one-to-one correspondence between spatial mapping and unit data.

[0093] Overlap analysis is performed on the geographic boundary contours of adjacent time phases. If the overlap is higher than the preset spatial continuity threshold, the geographic boundary contours of each time phase of the functional feature domain are integrated in time series to obtain the distribution of the functional feature domain in the stable time phase. If the overlap is lower than the spatial continuity threshold, the offset of the center position and the rate of change of the area of ​​the functional feature domain are calculated, and the offset of the center position and the rate of change of the area of ​​the functional feature domain are integrated in time series to obtain the distribution of the functional feature domain in the evolution phase.

[0094] In this embodiment, after obtaining the functional feature domain contours of consecutive time phases, an overlap analysis is performed on the contours of adjacent time phases. The overlap is obtained by calculating the ratio of the intersection area to the joint area of ​​the contour polygons, and a spatial consistency threshold is set. As a criterion for judgment. If the overlap is higher than... Then, adjacent contours are fused sequentially according to time, that is, linear interpolation and vertex integration are performed on the polygon nodes of each time phase to obtain a stable temporal functional feature domain contour. This contour represents the spatial distribution range of the functional feature domain that remains basically stable within the observation period, while retaining the unit number and cluster information of each unit to ensure data closure. If the overlap is lower than the threshold If this is the case, then the functional feature domain is considered to have undergone significant spatial evolution, requiring dynamic feature integration. First, the change in the centroid coordinates of this domain is calculated. and area change rate ,in It is used to quantify the offset of the center position and the trend of spatial expansion or contraction. , The functional feature domain is represented in the first... The coordinates of the centroid (center point) at each observation phase. Here Corresponding to east-west coordinates, Corresponding to north-south coordinates. The functional feature domain is represented in the first... The area under each observation phase is usually expressed in square meters. Area is used to quantify the spatial expansion or contraction of functional feature domains and is an important parameter for evolution trend analysis. Centroid calculation is typically obtained by averaging the coordinates of all nodes of the functional feature domain polygon, reflecting the spatial center position of the entire feature domain. Subsequent analysis of continuous time phases... and Perform weighted time series integration and set time weights. Correlation with the current time phase Make adjustments to give more weight to the most recent phase than to earlier phases (e.g.) , ), representing the time from the reference time (usually the current phase or the initial phase) to the 1st phase. The time intervals of each phase, which can be in hours or days, are used to obtain the contour of the functional characteristic domain of the evolution phase. This contour is also stored in a polygon node structure with accompanying evolution parameter annotations, which are used for subsequent ecological function score calculation and trend analysis to ensure that the temporal spatial changes are fully quantified.

[0095] Of particular importance is that the distribution of the functional characteristic domains of each time phase is the distribution of the functional characteristic domains of the stable time phase and the distribution of the functional characteristic domains of the evolution time phase.

[0096] Optionally, the evaluation of grassland ecosystem function scores in step S5 includes:

[0097] The grassland ecosystem function score of the target area is evaluated by using the area proportion and spatial continuity of the functional characteristic domains in the distribution of stable phases as the basic score, and the area change rate of the distribution of functional characteristic domains in the evolution phases and the center shift rate obtained from the center position shift as deduction factors.

[0098] In this embodiment, a basic score is calculated for each functional feature domain based on its distribution during stable time phases. The basic score comprehensively considers the area proportion and spatial continuity of the functional feature domain, and is obtained through a weighted average. The basic score for each functional feature domain can be described as: Basic Score = Area Proportion Contribution + Spatial Continuity Contribution. The area proportion is determined by the ratio of the functional feature domain area to the total area of ​​the target area, and spatial continuity is assessed by the degree of aggregation of vegetation distribution units within the functional feature domain. Subsequently, based on the distribution of functional feature domains during evolutionary time phases, a deduction factor is calculated, mainly including the area change rate and the center position shift rate. The area change rate reflects the magnitude of change in the functional feature domain area over time, and the center position shift rate reflects the movement of the centroid position over time. The deduction factor can be summarized as: Deduction Factor = Area Change Influence + Center Position Shift Influence. Finally, the basic score is subtracted from the deduction factor to obtain the comprehensive score for each functional feature domain, and then weighted by area proportion to obtain the overall grassland ecosystem function score of the target area.

[0099] Optionally, this specification also provides a grassland ecosystem function evaluation system based on multi-source data, used to perform the grassland ecosystem function evaluation method based on multi-source data as described above. The grassland ecosystem function evaluation system based on multi-source data includes:

[0100] The image registration module is used to acquire remote sensing image data and aerial images of the target area, and to perform geometric registration of the remote sensing image data and aerial images.

[0101] The spatial distribution extraction module is used to extract vegetation spatial distribution indicators of the target area based on the spectral characteristics of ground features from the geometric registration results.

[0102] The energy balance curve construction module is used to invert the surface temperature index based on the vegetation spatial distribution index, and combine the surface temperature index inversion results with the acquired meteorological observation data to construct the surface energy balance curve.

[0103] The functional feature domain segmentation module is used to identify the functional trends of each season based on the time series changes of the surface energy balance curve; to segment the functional feature domains according to the spatial clustering results of the functional trends of each season; and to determine the distribution of the functional feature domains in each time phase.

[0104] The functional score assessment module is used to evaluate the grassland ecosystem functional score based on the distribution of functional characteristic domains in each time phase.

[0105] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the application are intended to be included within the invention.

[0106] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A method for evaluating grassland ecosystem function based on multi-source data, characterized in that, Includes the following steps: Step S1: Acquire remote sensing image data and aerial images of the target area, and perform geometric registration of the remote sensing image data and aerial images; Step S2: Based on the spectral characteristics of ground features obtained from the geometric registration results, extract the spatial distribution indicators of vegetation in the target area; Step S3: Invert the surface temperature index based on the spatial distribution index of vegetation, and combine the surface temperature index inversion results with the acquired meteorological observation data to construct the surface energy balance curve; Step S3 involves constructing the surface energy balance curve, including: Acquire meteorological observation data; synchronize the surface temperature index with the meteorological observation data at the same timestamp, and pair them according to the spatial location corresponding to the vegetation distribution unit to obtain a multi-source data set of the unit; Input the multi-source data set of the unit into the preset energy budget equation to calculate the latent heat flux of each vegetation distribution unit; The temporal variation of latent heat flux in each vegetation distribution unit is plotted as a surface energy balance curve. Step S4: Based on the time series changes of the surface energy balance curve, identify the functional trends of each season; divide the functional characteristic domains according to the spatial clustering results of the functional trends of each season, and determine the distribution of the functional characteristic domains of each time phase. Step S5: Evaluate the grassland ecosystem function score based on the distribution of functional characteristic domains in each time phase.

2. The grassland ecosystem function evaluation method based on multi-source data according to claim 1, characterized in that, Step S1, which involves performing image geometric registration, includes: Detect the edge feature differences between any image frame in remote sensing image data and aerial images, and select control pixels based on the edge feature differences; Calculate the similarity matrix of control pixels in the image frame and the aerial image respectively to obtain the correspondence of control pixels; Based on the coordinate deviation in the correspondence of control pixels, the spatial mapping parameters between the image frame and the aerial image are determined respectively, and geometric resampling is performed on the pixel gray values ​​of the image frame and the aerial image according to the spatial mapping parameters. The geometric resampling results are integrated in time according to the timestamp, and the time integration result is used as the geometric registration result.

3. The grassland ecosystem function evaluation method based on multi-source data according to claim 2, characterized in that, Selecting control pixels includes: Edge detection and texture gradient calculation are performed on any image frame and aerial image in the remote sensing image data, respectively. The edge detection results and texture gradient information of the image frame and aerial image are then normalized and fused according to pixel position to obtain feature response pairs. Based on the texture repetition of each pixel in the feature response pair, regions in image frames and aerial images with texture repetition exceeding a preset texture repetition threshold are removed. The region removal results are divided into grids, and the cell with the highest response value in each grid is selected as a candidate control point. Calculate the local grayscale values ​​of candidate control points at different resolutions in remote sensing image data and aerial images, and select candidate control points whose local grayscale value standard deviation is less than the average grayscale value of the image as control pixels.

4. The grassland ecosystem function evaluation method based on multi-source data according to claim 1, characterized in that, Step S3 involves retrieving surface temperature indicators based on vegetation spatial distribution indices, including: By combining the infrared band information of remote sensing image data in the geometric registration results with vegetation spatial distribution indicators, the surface thermal radiation intensity of each vegetation distribution unit in the vegetation spatial distribution indicators is estimated. The spectral emissivity parameters of each vegetation distribution unit are calculated based on the surface thermal radiation intensity, and the surface thermal radiation intensity is corrected using the spectral emissivity parameters. The corrected surface thermal radiation intensity is converted into surface temperature, and a spatial weighted average is performed on the surface temperatures at the same time point to obtain the surface temperature index.

5. The grassland ecosystem function evaluation method based on multi-source data according to claim 1, characterized in that, Step S4 involves identifying the functional trends for each season, including: The spatial distribution indicators of vegetation and the surface energy balance curve at the same time stamp are spatiotemporally synchronized according to the vegetation distribution unit, and the spatiotemporal synchronization results are used as the data source to be analyzed. Based on the data sources to be analyzed, determine the dynamic trends of undercurrent heat flux and vegetation cover in each season; The dynamic trend of latent heat flux in each season is converted into the vegetation transpiration activity trend by the unit area of ​​each vegetation distribution unit. The vegetation transpiration activity trend and the vegetation cover dynamic trend in each season are then standardized and aggregated as the functional trend of each season.

6. The grassland ecosystem function evaluation method based on multi-source data according to claim 1, characterized in that, Step S4 involves dividing the functional feature domains, including: Based on the functional trends of each season, the trend of vegetation transpiration activity is taken as the first feature dimension, and the dynamic trend of vegetation coverage is taken as the second feature dimension. Calculate the trend similarity of the first feature dimension and the second feature dimension between each vegetation distribution unit, and perform hierarchical clustering based on the trend similarity, merging distribution unit clusters with similarity higher than the preset similarity threshold layer by layer; Based on the merging results of the distribution unit clusters in the first feature dimension and the second feature dimension, vegetation distribution units that simultaneously satisfy the conditions of spatial proximity and show a unidirectional change relationship in both the first and second feature dimensions are divided into functional feature domains.

7. The grassland ecosystem function evaluation method based on multi-source data according to claim 1, characterized in that, Step S4, determining the distribution of functional characteristic domains for each time phase, includes: The preset remote sensing observation period is used as the time phase, and the geographic boundary contour of each functional feature domain is extracted in the current time phase based on the spatial position relationship of vegetation distribution units in the functional feature domain. Overlap analysis is performed on the geographic boundary contours of adjacent time phases. If the overlap is higher than the preset spatial continuity threshold, the geographic boundary contours of each time phase of the functional feature domain are integrated in time series to obtain the distribution of the functional feature domain in the stable time phase. If the overlap is lower than the spatial continuity threshold, the offset of the center position and the rate of change of the area of ​​the functional feature domain are calculated, and the offset of the center position and the rate of change of the area of ​​the functional feature domain are integrated in time series to obtain the distribution of the functional feature domain in the evolution phase.

8. The grassland ecosystem function evaluation method based on multi-source data according to claim 1, characterized in that, Step S5, which evaluates grassland ecosystem function scores, includes: The grassland ecosystem function score of the target area is evaluated by using the area proportion and spatial continuity of the functional characteristic domains in the distribution of stable phases as the basic score, and the area change rate of the distribution of functional characteristic domains in the evolution phases and the center shift rate obtained from the center position shift as deduction factors.

9. A grassland ecosystem function evaluation system based on multi-source data, characterized in that, For performing the grassland ecosystem function assessment method based on multi-source data as described in claim 1, the grassland ecosystem function assessment system based on multi-source data includes: The image registration module is used to acquire remote sensing image data and aerial images of the target area, and to perform geometric registration of the remote sensing image data and aerial images. The spatial distribution extraction module is used to extract vegetation spatial distribution indicators of the target area based on the spectral characteristics of ground features from the geometric registration results. The energy balance curve construction module is used to invert the surface temperature index based on the vegetation spatial distribution index, and combine the surface temperature index inversion results with the acquired meteorological observation data to construct the surface energy balance curve. The functional feature domain segmentation module is used to identify the functional trends of each season based on the time series changes of the surface energy balance curve; to segment the functional feature domains according to the spatial clustering results of the functional trends of each season; and to determine the distribution of the functional feature domains in each time phase. The functional score assessment module is used to evaluate the grassland ecosystem functional score based on the distribution of functional characteristic domains in each time phase.

Images