Land ecological condition monitoring system and method based on remote sensing data

Abstract

The invention belongs to the field of remote sensing monitoring, and particularly relates to a land ecological condition monitoring system and method based on remote sensing data, and the method comprises the steps: obtaining a remote sensing image sequence through a preprocessing module, and obtaining initial cultivated land segmentation sub-regions and edge features; a feature extraction module extracts spatial features and fluctuation values such as normalized vegetation indexes and enhanced vegetation indexes based on the feature extraction module, and constructs a space-time coordinate system; the region determination module determines a real cultivated land area and a cultivated land area fluctuation value; the inversion module inverts the vegetation coverage, the moisture content and the soil fertility, and obtains an ecological degradation fluctuation value by combining with phenological period anomaly recognition; the early warning module performs real-time early warning based on the cultivated land area fluctuation value and the ecological degradation fluctuation value, and visually marks the position, the type and the severity of a change region through a GIS platform; according to the invention, precise monitoring and dynamic early warning of the ecological condition of the land are realized.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing monitoring, and in particular relates to a land ecological status monitoring system and method based on remote sensing data. Background Technology

[0002] As an important ecosystem, the assessment of the service value and ecological status of arable land is crucial for implementing arable land protection and natural resource management. Among existing technologies, Chinese patent CN118731935B discloses a method for calculating the ecosystem service value of arable land resources based on multi-source remote sensing technology. This method integrates temporal optical images, radar images, and field sample data to extract crop categories and planting areas, constructing a yield correction model coupled with topographic factors and a regional adaptive model of equivalent factors to calculate the ecosystem service value of arable land. Chinese patent application CN119761854A discloses a method for assessing the ecosystem services of arable land based on remote sensing monitoring. This method acquires multi-dimensional arable land data using remote sensing technology, builds an assessment model after dimensionality reduction, and achieves the assessment of arable land ecosystem services through dynamic weight allocation and service synergy optimization. Both approaches attempt to improve the scientific rigor of farmland ecological assessment through remote sensing technology, but they still face technical limitations in practical applications. Specifically, the farmland resource ecosystem service value accounting method based on multi-source remote sensing technology focuses on value accounting, but it lacks a real-time dynamic monitoring mechanism and fails to adequately integrate time-sensitive ecological parameters such as soil moisture and crop phenology, making it difficult to reflect instantaneous changes in farmland ecology. The farmland ecosystem service assessment method based on remote sensing monitoring focuses on ecosystem service assessment, relies on a single remote sensing data source, and fails to fully integrate multi-source heterogeneous data, such as radar and optical imagery. Furthermore, the assessment indicators are not linked to key ecological parameters such as soil fertility and surface temperature, resulting in incomplete assessment dimensions. Summary of the Invention

[0003] To address the shortcomings of existing technologies, this invention proposes a land ecological status monitoring system and method based on remote sensing data, comprising: a preprocessing module acquiring remote sensing image sequences and obtaining initial cultivated land sub-regions and edge features; a feature extraction module extracting spatial features and fluctuation values ​​such as normalized vegetation index and enhanced vegetation index, and constructing a spatiotemporal coordinate system; a region determination module determining the actual cultivated land area and cultivated land area fluctuation values; an inversion module inverting vegetation cover, moisture content, and soil fertility, and obtaining ecological degradation fluctuation values ​​by combining phenological anomaly identification; and an early warning module providing real-time early warning based on cultivated land area fluctuation values ​​and ecological degradation fluctuation values, and visually marking the location, type, and severity of changed areas through a GIS platform. This application achieves accurate monitoring and dynamic early warning of land ecological status.

[0004] To achieve the above objectives, the present invention provides the following technical solution:

[0005] The land ecological status monitoring system based on remote sensing data includes: a preprocessing module, a feature extraction module, a region determination module, an inversion module, and an early warning module;

[0006] The preprocessing module is used to acquire remote sensing image sequences of the target area and analyze them to obtain initial farmland segmentation sub-region sequences and edge segmentation feature sets;

[0007] The feature extraction module obtains the remote sensing image feature space based on the initial cultivated land segmentation sub-region sequence and edge segmentation feature set. The remote sensing image feature space includes the spatial feature sequence and spatial feature fluctuation value of any N edge extension directions of the remote sensing image. The spatial feature sequence of any N edge extension directions includes the normalized vegetation index sequence and enhanced vegetation index sequence of the corresponding extension direction. The spatial feature fluctuation value includes the spatial change value of the normalized vegetation index and the spatial change value of the enhanced vegetation index of the remote sensing image at each time point.

[0008] The region determination module determines the area of ​​the actual cultivated land sub-region based on the spatial features of remote sensing images and the initial cultivated land sub-regions, and determines the cultivated land area fluctuation value based on the determined changes in the actual cultivated land sub-region area.

[0009] The inversion module, based on the spatial features of remote sensing images, combines surface temperature with inversion to obtain the vegetation cover of cultivated land and obtain the ecological degradation fluctuation value of the real cultivated land sub-regions.

[0010] The early warning module provides real-time early warnings based on the fluctuation values ​​of cultivated land area or ecological degradation in the actual cultivated land sub-regions, combined with corresponding preset early warning values. It also automatically generates early warning information and visualizes and marks the spatial location, type, and degree of change in the actual cultivated land sub-regions through the GIS platform.

[0011] Specifically, the preprocessing module includes an acquisition unit and an image segmentation unit;

[0012] The acquisition unit is used to acquire a sequence of remote sensing images of the target area over a continuous preset time period, and to perform radiometric and geometric correction preprocessing on the remote sensing image sequence to obtain the corrected remote sensing image sequence.

[0013] The image segmentation unit is used to perform initial region identification based on the corrected remote sensing image sequence and the pre-trained YOLOv5 model. It obtains the initial region identification range and corresponding region type of remote sensing images acquired at different time points. Based on the initial region identification range and the image segmentation algorithm, it obtains the edge segmentation feature set corresponding to the initial cultivated land segmentation sub-region at each time point. The region types include cultivated land, construction land, forest land, pond, and abandoned land.

[0014] Specifically, the feature extraction module includes a center determination unit and a feature extraction unit;

[0015] The center determination unit obtains the first center point based on the pixel position points in the edge segmentation feature set corresponding to the initial cultivated land segmentation sub-region at each time point, combined with the minimum bounding rectangle, and obtains the centroid position as the second center point based on the edge pixel position points; and determines the true center position point of the initial cultivated land segmentation sub-region at each time point based on the average of the first center point and the second center point.

[0016] Two-dimensional coordinates are constructed based on the initial farmland sub-regions of the real center location at each time point, and the line connecting all aligned real center locations at consecutive time points is used as a three-dimensional time axis to construct a spatiotemporal coordinate system. The initial farmland sub-regions of the remote sensing images collected at each time point are embedded into the time point corresponding to the spatiotemporal coordinate system.

[0017] The feature extraction unit, based on the initial cultivated land sub-region at each time point, takes the true center point of the corresponding initial cultivated land sub-region as the starting point and performs feature extraction along any N edge extension directions in the initial cultivated land sub-region at the corresponding time point to obtain the spatial feature sequence and spatial feature fluctuation value of any N edge extension directions at the current time point.

[0018] Specifically, the region determination module includes a region determination unit, a region area unit, and a region fluctuation determination unit;

[0019] The region determination unit, based on the spatial feature sequence and spatial feature fluctuation value of any N edge extension directions at the current time point, combined with the region type fluctuation threshold, if the spatial feature fluctuation value of M consecutive pixel positions in any edge extension direction does not meet the region type fluctuation threshold, then the pixel position point corresponding to the largest spatial feature fluctuation value is taken as the edge point of the current edge extension direction, and this is repeated N times to determine the second edge position point set;

[0020] The regional area unit is based on the second edge position point set of the initial cultivated land sub-region at each time point and the corresponding initial cultivated land sub-region combined with the regional integration algorithm to obtain the area of ​​the actual cultivated land sub-region at each time point. At the same time, taking the pixel point corresponding to the largest actual cultivated land sub-region as the starting point, along the positive and negative axes of the three-dimensional time axis, the absolute value of the time feature change value of the area of ​​adjacent actual cultivated land sub-regions at the same pixel position point at consecutive time points is extracted.

[0021] Specifically, the region determination module includes a region fluctuation determination unit. The region fluctuation determination unit is used to obtain a first cultivated land area fluctuation value sequence based on the area of ​​the actual cultivated land sub-regions corresponding to adjacent time points, and at the same time, based on the absolute value of the time feature change value combined with the region type fluctuation threshold, obtain the proportion of the absolute value of the time feature change value of the same pixel location points that do not meet the region type fluctuation threshold, as the second cultivated land area fluctuation value sequence; and based on the average of the first cultivated land area fluctuation value and the second cultivated land area fluctuation value at corresponding adjacent time points in the first cultivated land area fluctuation value sequence and the second cultivated land area fluctuation value sequence, obtain the cultivated land area fluctuation value at adjacent time points.

[0022] Specifically, the inversion module includes a coverage inversion unit;

[0023] The vegetation cover inversion unit, based on the normalized vegetation index (NVI) sequence and enhanced vegetation index (EVI) sequence corresponding to each time point, and the spatial variation values ​​of the NVI and EVI at the corresponding time points, inverts to obtain the vegetation cover of each real cultivated land sub-region at the corresponding time point and the corresponding vegetation cover distribution probability function. At the same time, based on the vegetation cover of each real cultivated land sub-region and the absolute value of the temporal feature change value of the area of ​​adjacent real cultivated land sub-regions at the same pixel position at consecutive time points, the temporal distribution function of cultivated land vegetation cover is obtained.

[0024] Specifically, the inversion module also includes a first trend unit;

[0025] The first trend unit is used to invert the moisture content distribution state function of the corresponding real farmland sub-region based on the farmland vegetation coverage and the corresponding coverage distribution probability function of each real farmland sub-region, combined with the surface temperature and thermal inertia model. At the same time, based on the moisture content distribution state function of the real farmland sub-region combined with the farmland vegetation coverage time distribution function, the moisture content time fluctuation function of the real farmland sub-region at continuous time points is obtained.

[0026] Based on the vegetation cover and the corresponding cover distribution probability function of each real cultivated land sub-region, combined with the moisture content distribution state function of each real cultivated land sub-region, the soil fertility distribution function of each real cultivated land sub-region is obtained by inversion. At the same time, based on the soil fertility distribution function and the time distribution function of cultivated land vegetation cover of each real cultivated land sub-region at continuous time points, the time change trend of the first soil fertility at continuous time points is obtained by inversion.

[0027] Specifically, the inversion module also includes a second trend unit and a comprehensive trend unit;

[0028] The second trend unit, based on the spatial change values ​​of normalized vegetation index and enhanced vegetation index combined with the crop phenological extraction algorithm, obtains the crop phenological nodes corresponding to the real cultivated land segmented sub-regions. Based on the temporal feature change values ​​of the crop phenological nodes corresponding to the real cultivated land segmented sub-regions, it performs phenological anomaly identification, obtains the corresponding crop phenological anomaly identification feature space, and inverts the second soil fertility temporal change trend of the corresponding real cultivated land segmented sub-region based on the corresponding crop phenological anomaly identification feature space.

[0029] The integrated trend unit is used to time-align the first soil fertility time change trend with the second soil fertility time change trend, and then use a weighted average algorithm to combine the first soil fertility time change trend with the second soil fertility time change trend within the aligned time period to obtain the ecological degradation fluctuation value of the real cultivated land segmented sub-region.

[0030] Specifically, the early warning module includes an early warning discrimination unit and a real-time annotation unit;

[0031] The early warning judgment unit, based on the fluctuation values ​​of cultivated land area and ecological degradation at consecutive time points, determines that there is an anomaly in the sub-region of real cultivated land division at the corresponding consecutive adjacent time points if either the fluctuation value of cultivated land area or the fluctuation value of ecological degradation at two adjacent time points does not meet the corresponding early warning value.

[0032] The real-time annotation unit determines the spatial location and area of ​​the changed region based on the absolute value of the temporal characteristic change of adjacent real cultivated land sub-regions with anomalies at the same pixel location point under continuous time points, the fluctuation value of the first cultivated land area within the corresponding time period, and the corresponding pixel location point. Simultaneously, based on the normalized difference vegetation index (NDVI) and enhanced vegetation index (EGI) within the changed region and their corresponding spatial change values, it determines the type of change. Based on the determined spatial location and area of ​​the changed region, combined with the magnitude of the absolute value of the temporal characteristic change within the corresponding continuous time period, an evaluation algorithm is used to obtain the degree of change of cultivated land type within the current changed region. The spatial location, area of ​​change, type of change, and degree of change are then dynamically annotated and displayed in real-time on the real cultivated land sub-regions corresponding to adjacent time points in the spatiotemporal coordinate system using an automatic annotation algorithm.

[0033] Methods for monitoring land ecological conditions based on remote sensing data include:

[0034] Acquire remote sensing image sequences of the target area and analyze them to obtain initial farmland segmentation sub-region sequences and edge segmentation feature sets;

[0035] Based on the initial farmland segmentation sub-region sequence and edge segmentation feature set, the remote sensing image feature space is obtained;

[0036] Based on the spatial features of remote sensing images combined with the initial farmland sub-regions, the area of ​​the actual farmland sub-regions is determined, and the farmland area fluctuation value is determined based on the determined change value of the actual farmland sub-region area.

[0037] Based on the spatial characteristics of remote sensing images, combined with surface temperature and inversion, the vegetation cover of cultivated land is obtained, and the ecological degradation fluctuation value of the real cultivated land sub-regions is obtained.

[0038] Real-time early warnings are issued based on the fluctuation values ​​of cultivated land area or ecological degradation in real cultivated land sub-regions, combined with corresponding preset early warning values. At the same time, early warning information is automatically generated, and the spatial location, type, and degree of change of the real cultivated land sub-regions are visualized and marked through the GIS platform.

[0039] Compared with the prior art, the beneficial effects of the present invention are:

[0040] This invention addresses the shortcomings of existing technologies by constructing a multi-module collaborative remote sensing monitoring system. This system achieves high-precision dynamic perception and intelligent early warning of farmland ecological conditions. The preprocessing module improves the accuracy of initial farmland segmentation through radiometric and geometric correction combined with YOLOv5 region identification. The feature extraction module quantifies spatial feature fluctuations based on spatiotemporal coordinates and multi-directional vegetation index sequence analysis. The region determination module accurately captures farmland area changes through dynamic edge point calibration and area integration algorithms. The inversion module integrates vegetation cover, surface temperature, moisture content, and phenological characteristics to achieve multi-dimensional inversion of soil fertility and ecological degradation trends. The early warning module combines GIS visualization and automatic annotation technology to achieve real-time dynamic display of the location, type, and severity of changes in the area. This application effectively improves the automation, accuracy, and real-time performance of land ecological monitoring, providing reliable technical support for farmland protection and ecological governance. Attached Figure Description

[0041] Figure 1 This is a block diagram of the land ecological status monitoring system based on remote sensing data, as described in Embodiment 1 of the present invention.

[0042] Figure 2 This is a flowchart of the land ecological status monitoring method based on remote sensing data in Embodiment 2 of the present invention. Detailed Implementation

[0043] Example 1:

[0044] Please see Figure 1 The present invention provides an embodiment of a land ecological status monitoring system based on remote sensing data, comprising: a preprocessing module, a feature extraction module, a region determination module, an inversion module, and an early warning module;

[0045] The preprocessing module is used to acquire remote sensing image sequences of the target area and analyze them to obtain initial farmland segmentation sub-region sequences and edge segmentation feature sets. It should be further noted that the preprocessing module in this embodiment includes an acquisition unit and an image segmentation unit.

[0046] The acquisition unit is used to acquire a sequence of remote sensing images of the target area over a continuous preset time period, and to perform radiometric and geometric correction preprocessing on the remote sensing image sequence to obtain the corrected remote sensing image sequence.

[0047] It should be further explained that one specific implementation of obtaining the corrected remote sensing image sequence in this embodiment is as follows:

[0048] Based on the digital quantization value of remote sensing imagery, which is an integer with a preset number of bits and takes the value within a preset minimum to a preset maximum range, a sensor calibration algorithm is used to convert the digital quantization value into a radiance value received by the sensor. This radiance value takes the value within a preset radiance range and the accuracy meets the preset radiance accuracy standard.

[0049] Based on this radiance value, an atmospheric correction algorithm is used, combined with atmospheric parameters acquired during image acquisition. Specifically, the aerosol optical thickness is taken within a preset wavelength range, with accuracy meeting preset aerosol accuracy standards; the water vapor content is taken within a preset water vapor content range, with accuracy meeting preset water vapor accuracy standards. This eliminates the influence of atmospheric scattering and absorption on the radiance value, obtaining surface reflectance data. This surface reflectance data is taken within a preset reflectance range, with accuracy meeting preset reflectance accuracy standards, thus completing the radiometric correction.

[0050] Based on radiometrically corrected remote sensing images and high-precision reference maps, the reference maps use a preset coordinate system with coordinate accuracy conforming to preset map accuracy standards. Through a ground control point selection algorithm, no fewer than a preset number of evenly distributed ground control points with distinct characteristics are selected on the images and reference maps. The pixel coordinates of each ground control point in the images and the geographic coordinates in the reference maps are recorded, with pixel coordinate accuracy conforming to preset pixel accuracy standards and geographic coordinate accuracy conforming to preset geographic coordinate accuracy standards.

[0051] Based on the pixel coordinates and geographic coordinates of ground control points, a geometric transformation model algorithm is used to construct the transformation relationship between image pixel coordinates and reference map geographic coordinates using a quadratic polynomial. The coefficient parameters of the polynomial transformation model are calculated, and the accuracy of these coefficient parameters meets the preset coefficient accuracy standard.

[0052] Based on the coefficient parameters of the polynomial transformation model, the position of each pixel in the image is recalculated and rearranged using a bilinear interpolation method through a resampling algorithm. This ensures that the positional deviation of the resampled pixels does not exceed a preset pixel deviation threshold, thus ensuring that the coordinate system of the image is consistent with the coordinate system of the reference map and that the coordinate deviation does not exceed a preset coordinate deviation threshold, thereby obtaining a geometrically corrected remote sensing image sequence.

[0053] The image segmentation unit is used to perform initial region identification based on the corrected remote sensing image sequence and the pre-trained YOLOv5 model. It obtains the initial region identification range and corresponding region type of remote sensing images acquired at different time points. Based on the initial region identification range and the image segmentation algorithm, it obtains the edge segmentation feature set corresponding to the initial cultivated land segmentation sub-region at each time point. The region type includes, but is not limited to, cultivated land, construction land, forest land, pond, and abandoned land.

[0054] It should be further explained that one specific implementation of the image segmentation algorithm in this embodiment is as follows:

[0055] Based on the corrected remote sensing image sequence, a pre-trained YOLOv5 model was used. This model was trained on sample images of multiple regions, including cultivated land, construction land, forest land, ponds, and abandoned land. The corrected single-frame remote sensing image was input, and the feature extraction network of the model was used to perform multi-scale feature fusion on the image to obtain feature maps at different levels.

[0056] Based on the feature map, the YOLOv5 model's prediction head performs class probability prediction and bounding box regression on each preset anchor box, filters out anchor boxes with class probabilities greater than a preset probability threshold, merges anchor boxes with overlap exceeding a preset overlap threshold, obtains the initial recognition bounding box for each type of region, and determines the region type based on the highest class probability corresponding to the anchor box, thus completing the initial region recognition.

[0057] When the image segmentation unit performs image segmentation and edge feature extraction, based on the initial region recognition range and corresponding region type, it uses a semantic segmentation algorithm to refine the category labeling of image pixels within the initial region recognition range as the region of interest. Cultivated land is labeled as a set of pixels with continuous cultivation texture, construction land is labeled as a set of pixels with regular geometric shapes of artificial buildings, forest land is labeled as a set of pixels with dense vegetation texture, pond is labeled as a set of pixels with uniform water texture, and abandoned land is labeled as a set of pixels with no obvious cultivation traces and sparse vegetation cover.

[0058] Based on the refined annotation results, the gradient change at the boundary of pixels of different categories is calculated by the edge detection algorithm. Pixels with gradient values ​​exceeding the preset gradient threshold are extracted as edge pixels. Edge pixels within the same region type are connected according to spatial continuity to form a closed edge contour. The edge segmentation feature set corresponding to the initial cultivated land segmentation sub-region at each time point is obtained. This feature set contains the spatial coordinates of the edge pixels and the category difference information of adjacent pixels.

[0059] It should be further explained that in this embodiment, the gradient change at the boundary between pixels of different categories is calculated, and pixels with gradient values ​​exceeding a preset gradient threshold are extracted as edge pixels. The specific implementation method is as follows:

[0060] Based on the refined pixel category data, an edge detection algorithm is used to set a square neighborhood window of a preset size for each pixel in the image. This window contains the center pixel and a preset number of neighboring pixels. Based on the category information of each pixel in the window, the category of the center pixel is compared with the category of each neighboring pixel through a category difference calculation method. The same category is recorded as no difference, and different categories are recorded as different. The number of neighboring pixels with differences in the window is counted as the initial category difference value of the center pixel.

[0061] Based on the initial class difference value, the gradient components in the horizontal and vertical directions are calculated using the gradient operator: The horizontal gradient component is obtained by comparing the class difference between the center pixel and its left and right adjacent pixels. A positive value is assigned when the left neighbor pixel is different from the center pixel, and a negative value is assigned when the right neighbor pixel is different from the center pixel. The sum of the absolute values ​​of the two is taken as the horizontal gradient component. The vertical gradient component is obtained by comparing the class difference between the center pixel and its upper and lower adjacent pixels. A positive value is assigned when the upper neighbor pixel is different from the center pixel, and a negative value is assigned when the lower neighbor pixel is different from the center pixel. The sum of the absolute values ​​of the two is taken as the vertical gradient component.

[0062] Based on the gradient components in the horizontal and vertical directions, the gradient components in the two directions are combined and calculated using a gradient synthesis algorithm to obtain the gradient value of the center pixel. This gradient value changes positively with the intensity of the class difference.

[0063] The gradient value of each center pixel is compared with a preset gradient threshold. When the gradient value of the center pixel exceeds the preset gradient threshold, the center pixel is marked as a candidate edge pixel. Spatial continuity detection is performed on all candidate edge pixels. By judging whether the spatial distance between adjacent candidate edge pixels is within a preset range, candidate edge pixels with a distance within the preset range are grouped into the same edge segment. Then, the edge segments corresponding to the same region type are connected in spatial order to form a closed edge contour. The spatial coordinates of all pixels on the contour and the corresponding category difference values ​​are extracted as a component of the edge segmentation feature set.

[0064] The feature extraction module obtains the remote sensing image feature space based on the initial cultivated land segmentation sub-region sequence and edge segmentation feature set. The remote sensing image feature space includes the spatial feature sequence and spatial feature fluctuation value of any N edge extension directions of the remote sensing image. The spatial feature sequence of any N edge extension directions includes the normalized vegetation index sequence and enhanced vegetation index sequence of the corresponding extension direction. The spatial feature fluctuation value includes the spatial change value of the normalized vegetation index and the spatial change value of the enhanced vegetation index of the remote sensing image at each time point.

[0065] It should be further noted that the feature extraction module in this embodiment includes a center determination unit and a feature extraction unit;

[0066] It should be further explained that arable land often presents irregular shapes due to topography, farming methods, etc. Traditional center positioning methods are easily affected by shape, leading to deviations. Therefore, multiple algorithms are combined to dynamically adjust weights to accurately determine the true center, solving the problem of inaccurate center positioning of irregular arable land. At the same time, land ecological status monitoring needs to take into account both spatial characteristics and temporal dynamic changes. However, remote sensing image sub-regions at different time points lack a unified spatiotemporal reference. Therefore, a spatiotemporal coordinate system integrating spatial coordinates and temporal parameters is constructed, and sub-regions at each time point are embedded in it to achieve spatiotemporal correlation and unification.

[0067] The center determination unit obtains the first center point based on the pixel position points in the edge segmentation feature set corresponding to the initial cultivated land segmentation sub-region at each time point, combined with the minimum bounding rectangle, and obtains the centroid position as the second center point based on the edge pixel position points; and determines the true center position point of the initial cultivated land segmentation sub-region at each time point based on the average of the first center point and the second center point.

[0068] It should be further explained that one specific method for determining the true center location point in this embodiment is as follows:

[0069] Based on the spatial coordinates of edge pixels in the edge segmentation feature set corresponding to the initial cultivated land segmentation sub-region at each time point, the minimum area bounding rectangle algorithm is used to perform multi-directional rotation fitting on all edge pixels. It traverses multiple rotation directions within a preset angle range, calculates the area of ​​the rectangle that can contain all edge pixels in each direction, and selects the rectangle with the smallest area as the target bounding rectangle. The side of this rectangle can form a preset angle with the horizontal direction. The coordinates of the four vertices of the target bounding rectangle are extracted, and the horizontal and vertical coordinates of the intersection of the two diagonals are taken by the diagonal midpoint calculation method to form the coordinates of the first center point, so that the center point is more in line with the geometric center of the irregular cultivated land.

[0070] Based on the spatial coordinates of edge pixels and the class difference information of adjacent pixels in the edge segmentation feature set, a weighted centroid algorithm is used to assign weight values ​​to edge pixels. Edge pixels that are continuous with surrounding pixels and have stable class differences are given higher weights, while isolated edge pixels or edge pixels with abrupt class differences are given lower weights. The horizontal coordinates of all edge pixels are multiplied by their corresponding weights, summed, and then divided by the sum of all weights to obtain the weighted average horizontal coordinates. The vertical coordinates of all edge pixels are multiplied by their corresponding weights, summed, and then divided by the sum of all weights to obtain the weighted average vertical coordinates. The two are combined to form the coordinates of the second centroid, reducing the interference of abnormal pixels in irregular edges on the centroid.

[0071] Based on the coordinates of the first center point, the coordinates of the second center point, and the distribution characteristics of edge pixels, a dynamic weighting algorithm is used to calculate the average distance from the edge pixels to each side of the minimum area bounding rectangle. The larger the distance, the more irregular the shape of the farmland. In this case, the weight of the second center point is increased, and the weight of the first center point is increased if the distance is smaller. According to the dynamically determined weights, the horizontal coordinates of the first and second center points are summed according to their weights to obtain the horizontal coordinates of the true center. The vertical coordinates of the two are summed according to their weights to obtain the vertical coordinates of the true center. These are combined to form the true center location of the irregular initial farmland sub-region at each time point.

[0072] By employing a minimum area bounding rectangle algorithm to perform multi-directional rotation fitting on the edge pixels of irregular farmland, and selecting the midpoint of the diagonal of the minimum area bounding rectangle as the first center point, the system effectively avoids misjudging the geometric center of irregular areas, making the first center point more closely match the actual geometric distribution of the farmland. A weighted centroid algorithm is used to assign differentiated weights to edge pixels, reducing the interference from isolated or abruptly changing edge pixels and improving the representativeness of the second center point coordinates for edge distribution features. Combined with a dynamic weighting algorithm, the weight ratio of the two algorithms is adaptively adjusted according to the degree of irregularity in the farmland shape, ensuring that the final determined true center location simultaneously considers both geometric center and edge distribution features. This significantly improves the accuracy and stability of center positioning for irregular farmland, providing a reliable benchmark for subsequent features extraction based on center location and spatiotemporal coordinate system construction, thereby enhancing the accuracy and reliability of the entire land ecological status monitoring system.

[0073] Two-dimensional coordinates are constructed based on the initial farmland sub-regions of the real center location at each time point, and the line connecting all aligned real center locations at consecutive time points is used as a three-dimensional time axis to construct a spatiotemporal coordinate system. The initial farmland sub-regions of the remote sensing images collected at each time point are embedded into the time point corresponding to the spatiotemporal coordinate system.

[0074] It should be further explained that one specific method for constructing the spatiotemporal coordinate system in this embodiment is as follows:

[0075] Based on the true center location of the initial cultivated land sub-region at each time point, with the true center location as the origin of the coordinate system, the direction consistent with the horizontal scanning direction of the remote sensing image is set as the positive direction of the first coordinate axis, and the direction perpendicular to the horizontal scanning direction is set as the positive direction of the second coordinate axis, and a local two-dimensional coordinate system for the corresponding time point is established through the coordinate system orientation rules.

[0076] Based on this local two-dimensional coordinate system, by subtracting the coordinates of the real center point from the spatial coordinates of all edge pixels in the initial cultivated land sub-region using the relative coordinate calculation method, the relative coordinate values ​​of each edge pixel in the local two-dimensional coordinate system are obtained, forming a set of relative coordinates of edge pixels;

[0077] Based on the true center location points and corresponding time information at continuous time points, all true center location points are sorted according to the chronological order of remote sensing image acquisition using a time series sorting algorithm to form an ordered true center sequence.

[0078] Using a spatial connection algorithm, the true center locations of two adjacent time points in an ordered true center sequence are connected by straight line segments to form a continuous polyline as the spatial path of the three-dimensional time axis;

[0079] By using the time axis assignment rules, the acquisition time of each time point is used as the third axis parameter and is evenly assigned to the corresponding position on the three-dimensional time axis according to the time interval, so that each node of the three-dimensional time axis contains both spatial coordinates and time parameters.

[0080] Based on the local two-dimensional coordinate system at each time point, the spatial coordinates and time parameters of the corresponding nodes on the three-dimensional time axis, the origin of the local two-dimensional coordinate system is accurately mapped to the node position of the corresponding time parameter on the three-dimensional time axis through a coordinate mapping algorithm; through a plane adaptation algorithm, the plane containing the first and second coordinate axes of the local two-dimensional coordinate system is kept perpendicular to the tangent direction of the three-dimensional time axis at that node, ensuring a continuous transition of the local plane direction at different time points.

[0081] Based on the relative coordinate set of edge pixels, the relative coordinates of each edge pixel are superimposed with the spatial coordinates of the corresponding node on the three-dimensional time axis through an absolute coordinate transformation algorithm to obtain the absolute coordinates of each edge pixel in the three-dimensional spatiotemporal coordinate system.

[0082] Based on the absolute coordinates of all edge pixels, the complete outline of the initial farmland sub-region is restored in the vertical plane of the corresponding time point in the three-dimensional spatiotemporal coordinate system through the region reconstruction algorithm, thus completing the embedding of the sub-region at that time point. The above mapping, adaptation, transformation and reconstruction process is repeated for the initial farmland sub-regions at consecutive time points, and finally an integrated spatiotemporal coordinate system in which each sub-region is distributed in chronological order is formed.

[0083] This process effectively solves the technical challenge of inaccurate geometric center positioning in irregular farmland by integrating the minimum area bounding rectangle algorithm and the weighted centroid algorithm, and introducing a dynamic weight adjustment mechanism. The minimum area bounding rectangle algorithm accurately captures the macroscopic geometric shape of farmland through multi-directional rotation fitting, ensuring that the first center point conforms to the actual distribution. The weighted centroid algorithm suppresses abnormal edge pixel interference through differentiated weights, improving the ability of the second center point's coordinates to represent edge details. Furthermore, the weights of both algorithms are adaptively allocated according to the degree of shape irregularity, so that the final determined true center location point has both geometric inclusiveness and distribution representativeness, significantly improving the accuracy and robustness of center positioning. On this basis, by constructing a three-dimensional spatiotemporal coordinate system that integrates spatial coordinates and time dimensions, multi-temporal farmland sub-regions are uniformly embedded into the spatiotemporal framework, realizing the temporal serialization and spatial integrated management of monitoring data. This spatiotemporal coordinate system, through the combination of local two-dimensional coordinates with a three-dimensional time axis and planar vertical adaptation, ensures the spatial comparability and temporal continuity of multi-temporal data, providing a stable and reliable spatiotemporal benchmark for subsequent analyses such as feature extraction and change detection. Ultimately, this technical solution enhanced the system's ability to perceive and monitor changes in irregular farmland by improving the accuracy of central positioning and establishing a unified spatiotemporal benchmark.

[0084] The feature extraction unit, based on the initial cultivated land sub-region at each time point, takes the true center point of the corresponding initial cultivated land sub-region as the starting point and performs feature extraction along any N edge extension directions in the initial cultivated land sub-region at the corresponding time point to obtain the spatial feature sequence and spatial feature fluctuation value of any N edge extension directions at the current time point.

[0085] It should be further explained that, in this embodiment, one method for feature extraction along any N edge extension directions of the initial cultivated land segmentation sub-region at the corresponding time point is as follows:

[0086] Based on the edge segmentation feature set of the initial cultivated land sub-region at each time point, the edge segmentation feature set contains the spatial coordinates of edge pixels and the class difference information of adjacent pixels. The gradient direction of each edge pixel is quantized by the directional clustering algorithm, where the gradient direction is the direction from the current pixel to the adjacent but different pixel.

[0087] Statistically analyze the gradient direction distribution of all edge pixels, group edge pixels with direction differences within a preset range into the same direction cluster, and calculate the number of edge pixels contained in each direction cluster;

[0088] Select the direction cluster with the largest number of pixels, and take the average value of the gradient direction of each direction cluster as the representative direction of the cluster; perform angle calibration on the representative direction so that the angle difference between any two adjacent representative directions does not exceed the preset angle deviation, and obtain the angle parameters of any N edge extension directions.

[0089] Based on the initial farmland sub-region and the corresponding real center point at each time point, a ray is generated along each determined edge extension direction using the real center point as the ray starting point and the ray generation algorithm. The extension length of the ray is limited by the maximum radius of the initial farmland sub-region to ensure that the ray penetrates from the center to the edge.

[0090] The pixel traversal algorithm detects pixels sequentially along the ray direction, recording the spatial coordinates and sub-region attributes of each pixel traversed by the ray. Pixels that the ray touches after passing through the initial farmland sub-region are removed, forming a pixel sequence that is continuously arranged from the center to the edge in each direction. The sequence length is the total number of pixels traversed by the ray.

[0091] Based on the pixel sequence in each direction and the remote sensing image at the corresponding time point, the spectral reflectance of each pixel in the pixel sequence is processed by an index calculation algorithm. The normalized vegetation index is obtained by calculating the reflectance of the near-infrared band and the red band, while the enhanced vegetation index is obtained by calculating the reflectance of the near-infrared band, the red band, and the blue band. The normalized vegetation index value and the enhanced vegetation index value of each pixel are recorded sequentially according to the arrangement order of the pixels in the sequence, forming a normalized vegetation index sequence and an enhanced vegetation index sequence for each edge extension direction. The numerical position in the sequence corresponds one-to-one with the position of the corresponding pixel in the ray direction.

[0092] Based on the normalized vegetation index (NVI) sequence and the enhanced vegetation index sequence in each direction, the NVI values ​​of two adjacent pixels in the sequence are differentially calculated using a difference calculation algorithm. The absolute value of the result is taken as the spatial change value of the NVI of the adjacent pixel pair.

[0093] Perform the same difference operation on the enhanced vegetation index values ​​of two adjacent pixels in the sequence, and take the absolute value of the operation result as the spatial change value of the enhanced vegetation index of the adjacent pixel pair.

[0094] Arrange all spatial variation values ​​sequentially according to the order in which adjacent pixel pairs appear in the pixel sequence to form a spatial feature fluctuation value sequence that matches the length of the spatial feature sequence. Each fluctuation value corresponds to the difference between two adjacent pixels in the spatial feature sequence. Finally, obtain the spatial feature sequence and spatial feature fluctuation value contained in the feature space of the remote sensing image at the current time point.

[0095] This process significantly improves the comprehensiveness and accuracy of capturing spatial features of cultivated land vegetation by combining multi-directional ray scanning with gradient directional clustering. Edge pixel gradient directional clustering analysis autonomously determines the most representative edge extension direction, overcoming the subjectivity and limitations of manual direction selection and ensuring a high degree of consistency between the extracted feature direction and the actual contour features of the cultivated land. By emitting rays from the actual center point to each edge direction and traversing the pixel sequence, gradual feature acquisition from the core to the edge is achieved, effectively capturing the spatial heterogeneity and transitional features of cultivated land vegetation. The normalized vegetation index sequence and enhanced vegetation index sequence calculated based on the ray pixel sequence not only reflect the spatial distribution patterns of vegetation density and activity but also accurately quantify the boundary abrupt changes and internal gradual transition information of vegetation features through spatial feature fluctuation values ​​generated by differences between adjacent pixels. This multi-directional, multi-indicator, and multi-scale feature extraction system can comprehensively characterize the spatial distribution patterns and changing trends of cultivated land vegetation, providing high-precision and high-reliability feature data support for subsequent monitoring of cultivated land area changes and ecological degradation inversion, and greatly enhancing the system's perception and early warning accuracy for complex-shaped cultivated land and gradual ecological changes.

[0096] The region determination module determines the area of ​​the actual cultivated land sub-region based on the spatial features of remote sensing images and the initial cultivated land sub-regions, and determines the cultivated land area fluctuation value based on the determined changes in the actual cultivated land sub-region area. It should be further noted that the region determination module in this embodiment includes a region determination unit, a region area unit, and a region fluctuation determination unit.

[0097] The region determination unit, based on the spatial feature sequence and spatial feature fluctuation value of any N edge extension directions at the current time point, combined with the region type fluctuation threshold, if the spatial feature fluctuation value of M consecutive pixel positions in any edge extension direction does not meet the region type fluctuation threshold, then the pixel position corresponding to the largest spatial feature fluctuation value is taken as the edge point of the current edge extension direction. Based on this, it is repeated N times to determine the second edge position point set.

[0098] It should be further explained that the specific process of determining the second edge location point set in this embodiment includes:

[0099] Based on historical remote sensing images of the same growing season over several consecutive years within the target area, the images are divided into sub-regions such as cultivated land, construction land, forest land, ponds, and abandoned land using a regional type classification algorithm. Spatial feature fluctuation values ​​of all edge extension directions within each sub-region are extracted to form a categorized historical sample set.

[0100] For a sample set of the same region type, the arithmetic mean of the spatial characteristic fluctuation values ​​is calculated using a statistical analysis algorithm. At the same time, the degree of deviation of all values ​​in the sample set from the mean is calculated to obtain a statistical quantity reflecting the discrete characteristics.

[0101] Based on the mean and discrete statistics, combined with the inherent attributes of the region type, a threshold calibration algorithm is used to form the spatial characteristic fluctuation threshold for the region type. The lower limit is the discrete statistics obtained by subtracting a preset multiple from the mean, and the upper limit is the discrete statistics obtained by adding a preset multiple to the mean. The lower limit of the fluctuation threshold for cultivated land must be higher than the typical fluctuation value of abandoned land, and the upper limit must be lower than the abrupt fluctuation value of construction land. For example, the fluctuation threshold of the inherent attributes of the region type, such as the seasonal fluctuation of cultivated land due to crop growth, must cover the normal range of change from the greening period to the maturity period. Construction land, due to its stable structure, has a smaller fluctuation threshold set to limit abnormal changes.

[0102] Based on the spatial feature fluctuation value sequence of any N edge extension directions at the current time point (arranged in order from the real center point to the edge) and the corresponding fluctuation threshold of the cultivated land type, the algorithm checks the sequence sequentially starting from the first pixel position point through a point-by-point traversal. When the spatial feature fluctuation value of a certain pixel is less than the lower limit of the fluctuation threshold or greater than the upper limit of the fluctuation threshold, it is marked as a pixel that does not meet the fluctuation threshold, and continuous counting is started at the same time.

[0103] If subsequent adjacent pixels are consecutive pixels that do not meet the fluctuation threshold, the count is incremented until the consecutive count reaches the preset M pixels, and the spatial feature fluctuation values ​​of the M pixels do not fall within the fluctuation threshold range, then the traversal stops, the M pixel positions are arranged into a continuous subsequence, and the start and end positions of the subsequence are recorded.

[0104] Based on the spatial feature fluctuation values ​​of M consecutive locked pixel locations, a point-by-point comparison algorithm is used to compare the fluctuation value of each pixel in the subsequence in turn, and the pixel location with the largest fluctuation value is selected.

[0105] If multiple pixels have the same fluctuation value and all of them are the maximum value, then the pixel position closest to the edge direction in the subsequence is selected as the target point; this target point is marked as the edge point in the current edge extension direction, because its maximum fluctuation value corresponds to the critical abrupt change position of the internal features of the region to the external features, which can accurately define the boundary of the cultivated land sub-region.

[0106] The above detection, locking, and filtering operations are performed on any N edge extension directions, and an edge point is obtained in each extension direction. The spatial coordinate verification algorithm is used to check whether the spatial distribution of all edge points conforms to the contour trend of the initial cultivated land sub-region. Abnormal edge points that deviate from the overall trend are removed. If there are no edge points that meet the conditions in a certain direction, they are supplemented by interpolation based on the coordinates of edge points in adjacent directions. All the verified edge points are connected in sequence according to the angular order of the edge extension direction to form a closed edge contour. The spatial coordinates of all edge points contained therein are summarized into the second edge position point set.

[0107] This process effectively improves the accuracy and adaptability of regional edge detection by comprehensively utilizing historical remote sensing image data and spatial feature fluctuation analysis. Its core lies in constructing fluctuation threshold ranges for regional types based on multi-temporal images, fully considering the changing characteristics of different land cover attributes, ensuring that edge judgment is both statistically significant and conforms to the evolutionary patterns of land cover. By detecting consecutive pixel sequences that do not meet the fluctuation threshold in each direction and locating fluctuation extreme points, it can accurately capture the locations of abrupt feature changes, significantly suppressing misjudgments caused by noise or local variations, and ensuring that the edge points are located in the true transition zone between internal features and the external environment. Furthermore, by employing spatial verification and interpolation mechanisms, the spatial consistency and contour closure of edge points are guaranteed. The resulting second edge location point set not only has clear boundaries and a complete structure but also possesses high land use classification capabilities, making it particularly suitable for the precise delineation of easily confused areas such as cultivated land, abandoned land, and construction land. This provides reliable technical support for remote sensing image segmentation and dynamic land use monitoring.

[0108] The regional area unit is based on the second edge position point set of the initial cultivated land sub-region at each time point and the corresponding initial cultivated land sub-region combined with the regional integration algorithm to obtain the area of ​​the actual cultivated land sub-region at each time point. At the same time, taking the pixel point corresponding to the largest actual cultivated land sub-region as the starting point, along the positive and negative axes of the three-dimensional time axis, the absolute value of the time feature change value of the area of ​​adjacent actual cultivated land sub-regions at the same pixel position point at consecutive time points is extracted.

[0109] It should be further explained that the specific process of obtaining the area of ​​the actual cultivated land sub-region in this embodiment includes:

[0110] Based on the second edge location point set of the initial cultivated land segmentation sub-region at each time point and the corresponding edge segmentation feature set of the initial cultivated land segmentation sub-region, a boundary verification algorithm is used to compare the spatial position of each edge point in the second edge location point set with the edge pixels of the initial cultivated land segmentation sub-region. Abnormal edge points that exceed the range of the initial sub-region are removed, and valid edge points within the range of the initial sub-region are retained. The valid edge points are smoothed according to their spatial distribution to eliminate abrupt jumps between adjacent points and form a continuous set of boundary transition points.

[0111] Based on the smoothed effective edge points, a polygon fitting algorithm is used to connect all effective edge points in a clockwise or counterclockwise order according to the edge extension direction to form a closed boundary contour that surrounds the farmland segmentation sub-region, ensuring that the contour has no intersection and completely covers the internal farmland pixels.

[0112] Based on the closed boundary contour and the pixel resolution of the remote sensing image, the region within the closed boundary contour is divided into grid units with the same pixel size as the image through a grid division algorithm. The actual area of ​​each grid unit is determined by the pixel resolution. For example, if the pixel resolution is a preset length, then the area of ​​a single grid unit is the square of that length.

[0113] When performing region integration, the region integration algorithm is used to traverse and count all grid cells within the closed boundary contour. The number of grid cells completely contained within the contour is counted. For grid cells partially contained at the edge of the contour, the area proportion within the contour is calculated. For example, if the proportion exceeds a preset threshold, it is counted as a complete cell; otherwise, it is counted as a partial cell.

[0114] The area of ​​all fully contained grid cells is summed with the area of ​​the converted partial cells to obtain the actual farmland sub-region area corresponding to each time point.

[0115] It should be further explained that, in this embodiment, the absolute value of the temporal feature change of adjacent real farmland sub-regions at the same pixel location point along the positive and negative axes of the three-dimensional time axis is extracted. The specific process includes:

[0116] Based on the area of ​​the real farmland sub-regions at continuous time points, the extreme value screening algorithm is used to compare the area size at all time points, and the real farmland sub-region with the largest area is selected as the benchmark sub-region. The spatial coordinates of all pixel locations within the benchmark sub-region are recorded to form a benchmark pixel set.

[0117] Based on the time parameter sorting of the three-dimensional time axis, the time point corresponding to the reference sub-region is taken as the origin. Through the direction definition rules, the direction extending from the origin to the past time point is set as the negative axis direction of the three-dimensional time axis, and the direction extending to the future time point is set as the positive axis direction, thus clarifying the arrangement order of continuous time points on the positive and negative axes.

[0118] Based on each pixel location in the baseline pixel set, a coordinate mapping algorithm is used to map its spatial coordinates to the actual farmland sub-regions corresponding to adjacent time points along the positive and negative axes of the three-dimensional time axis. Pixel locations in each adjacent sub-region that have the same spatial coordinates as those in the baseline pixel set are found, forming a set of correspondences of the same pixel location at different time points. Pixels that cannot be matched with the same coordinates are marked as missing pixels.

[0119] Based on the positive and negative axes of the three-dimensional time axis, a sequence traversal algorithm is used to select the previous adjacent time points along the negative axis and the subsequent adjacent time points along the positive axis, starting from the time point corresponding to the reference sub-region, forming several pairs of adjacent time points, including the reference time point and the previous time point, the reference time point and the subsequent time point, the adjacent previous time points, and the adjacent subsequent time points.

[0120] Based on the correspondence of the same pixel position in adjacent time point pairs, the area ratio of the pixel position in the real farmland sub-region corresponding to the two adjacent time points is extracted by the difference operation algorithm. That is, the area ratio contributed by the pixel in its sub-region. The area ratio difference is obtained by subtracting the area ratio of the previous time point from the area ratio of the later time point.

[0121] Based on the absolute value conversion algorithm, the above area ratio difference is processed, and its absolute value is taken as the absolute value of the time feature change value of the pixel position point in the corresponding adjacent time point pair; for the position marked as missing pixel point, a preset missing value is assigned as the absolute value of the time feature change value, and finally a set of absolute values ​​of the time feature change values ​​of all the same pixel position points under continuous time points is formed.

[0122] The regional fluctuation determination unit is used to obtain a first cultivated land area fluctuation value sequence based on the area of ​​the actual cultivated land sub-regions corresponding to adjacent time points. At the same time, based on the absolute value of the time feature change value and the regional type fluctuation threshold, it obtains the proportion of the absolute value of the time feature change value of the same pixel location points that do not meet the regional type fluctuation threshold, as the second cultivated land area fluctuation value sequence. Based on the average of the first cultivated land area fluctuation value and the second cultivated land area fluctuation value at corresponding adjacent time points in the first cultivated land area fluctuation value sequence and the second cultivated land area fluctuation value sequence, it obtains the cultivated land area fluctuation value at adjacent time points.

[0123] It should be further explained that the specific process for obtaining the fluctuation value of cultivated land area at adjacent time points in this embodiment includes:

[0124] Based on the actual cultivated land sub-region area corresponding to continuous time points and adjacent time point pairs (each time point is paired with the next time point in chronological order on the three-dimensional time axis), the difference calculation method is used to calculate the difference between the actual cultivated land sub-region area of ​​the next time point and the actual cultivated land sub-region area of ​​the previous time point in each pair of adjacent time points. The absolute value of the difference is taken as the first cultivated land area fluctuation value of the adjacent time point pair. According to the arrangement order of the adjacent time point pairs on the three-dimensional time axis (from the earliest time point pair to the latest time point pair), all the first cultivated land area fluctuation values ​​are arranged in sequence to form the first cultivated land area fluctuation value sequence. Each fluctuation value in the sequence corresponds one-to-one with the corresponding adjacent time point pair.

[0125] Based on the set of absolute values ​​of temporal feature changes at the same pixel location point at consecutive time points (this set is extracted from the positive and negative axes of the three-dimensional time axis by the area unit) and the regional type fluctuation threshold corresponding to the cultivated land type (this threshold is obtained by statistical calibration of the spatial feature fluctuation values ​​of historical cultivated land samples), a threshold comparison algorithm is used to check the absolute values ​​of temporal feature changes of all pixel locations in each adjacent time point pair. If the absolute value of the temporal feature change of a pixel location point is greater than the upper limit of the regional type fluctuation threshold or less than its lower limit, then the pixel location point is determined to not meet the regional type fluctuation threshold. The number of pixel locations that do not meet the regional type fluctuation threshold in the adjacent time point pair is counted using a proportion statistics algorithm. The ratio of this number to the total number of pixel locations in the adjacent time point pair is taken as the second cultivated land area fluctuation value of the adjacent time point pair. The total number of pixel locations is the matching point of all pixels in the base sub-region in the adjacent sub-region, excluding missing pixels. All the second cultivated land area fluctuation values ​​are arranged in sequence according to the arrangement order of the adjacent time point pairs on the three-dimensional time axis to form a second cultivated land area fluctuation value sequence. Each fluctuation value in the sequence corresponds one-to-one with the corresponding adjacent time point pair.

[0126] Based on the first and second cultivated land area fluctuation value sequences (both sequences have the same length and the same number of adjacent time point pairs), the average value calculation method is used to sum the first and second cultivated land area fluctuation values ​​corresponding to the same adjacent time point pair in the two sequences. The summation result is then divided by two to obtain the average fluctuation value of the adjacent time point pair. This average fluctuation value is used as the cultivated land area fluctuation value of the adjacent time point. The above average calculation process is repeated for all adjacent time point pairs to obtain a set of cultivated land area fluctuation values ​​of consecutive adjacent time points. Each fluctuation value in the set corresponds to a set of adjacent time point pairs on the three-dimensional time axis.

[0127] This process effectively achieves high-precision dynamic monitoring of cultivated land area and its fluctuations by integrating multi-temporal spatial boundary information and time-series change characteristics. Its advantage lies in utilizing a dual verification mechanism of real edge point sets and the initial region, significantly improving the accuracy of area calculation and ensuring that the integration results more closely match the actual distribution of land features. By introducing pixel-level change tracking in the positive and negative directions of the three-dimensional time axis, it can keenly capture the differences in the area proportion of the same location at different time points, thereby extracting physically meaningful temporal feature change values. Furthermore, by combining macroscopic area fluctuations with microscopic pixel changes, a comprehensive fluctuation index is constructed through mean fusion. This index reflects both the overall area change of the region and incorporates information on abnormal changes in internal pixels, making the final cultivated land area fluctuation value both robust and sensitive. This method effectively suppresses misjudgments caused by transient interference or local noise, accurately identifying continuous changes caused by actual farming activities, and providing reliable data for dynamic monitoring of cultivated land resources, early warning of abandoned land, and land use assessment.

[0128] The inversion module, based on the spatial features of remote sensing images, combines surface temperature with inversion to obtain the vegetation cover of cultivated land and obtain the ecological degradation fluctuation value of the real cultivated land sub-regions.

[0129] It should be further noted that the inversion module in this embodiment includes a coverage inversion unit, a first trend unit, a second trend unit, and a comprehensive trend unit;

[0130] The vegetation cover inversion unit, based on the normalized vegetation index (NVI) sequence and enhanced vegetation index (EVI) sequence corresponding to each time point, and the spatial variation values ​​of the NVI and EVI at the corresponding time points, inverts to obtain the vegetation cover of each real cultivated land sub-region at the corresponding time point and the corresponding vegetation cover distribution probability function. At the same time, based on the vegetation cover of each real cultivated land sub-region and the absolute value of the temporal feature change value of the area of ​​adjacent real cultivated land sub-regions at the same pixel position at consecutive time points, the temporal distribution function of cultivated land vegetation cover is obtained.

[0131] It should be further explained that the inversion implementation process corresponding to the coverage inversion unit in this embodiment includes:

[0132] Based on a sample area with known arable land vegetation cover, where the sample area must include arable land plots with high, medium and low cover levels and have field-measured arable land vegetation cover data, the normalized vegetation index value, enhanced vegetation index value and corresponding measured cover value of each pixel in the sample area are extracted through a sample calibration algorithm to form an index-coverage sample set containing three types of parameters.

[0133] Based on the index-coverage sample set, the measured coverage values ​​in the sample set are divided into several continuous and non-overlapping coverage intervals according to the interval rule. The interval value is determined according to the range of sample coverage and the number of samples to ensure that the number of samples in each layer is balanced. Each interval constitutes a coverage layer, which contains the normalized vegetation index value, enhanced vegetation index value and measured coverage value of all pixels in the interval.

[0134] When assigning weights to samples in each layer, a weight allocation algorithm is used to assign a base weight to the normalized vegetation index (NDI) values ​​of all pixels within the layer: when the coverage value within the layer is in the low to medium range, the base weight is higher, reflecting the dominant influence of vegetation density on low-coverage areas; when the coverage value within the layer is in the high range, the base weight is lower, because the activity effect is more significant in high-density vegetation. A correction weight is assigned to the enhanced vegetation index values: the correction weight is complementary to the base weight, meaning the sum of the two is a fixed value. Layers with high coverage have higher correction weights, and layers with low coverage have lower correction weights, to reflect the corrective effect of vegetation activity on high-coverage areas after removing soil and atmospheric interference.

[0135] When fitting the local mapping relationship of each layer using the least squares method, the normalized vegetation index value of the pixels in the layer is multiplied by the base weight, and the enhanced vegetation index value is multiplied by the correction weight to obtain two weighted index values. Using these two weighted index values ​​as independent variables and the measured coverage value of the corresponding pixels as dependent variables, the coefficients of the linear regression equation are solved by the least squares method to obtain the local mapping relationship of each layer. This relationship is expressed in the form of an equation: predicted vegetation coverage value = a × weighted normalized vegetation index value + b × weighted enhanced vegetation index value + constant term, where a and b are coefficients, and a, b and the constant term are determined by fitting the samples in the layer.

[0136] When forming a global nonlinear mapping model based on the relationships between layers, an interlayer transition algorithm is used to select overlapping samples in the boundary interval of adjacent coverage layers and calculate the prediction deviation of the local mapping relationship between adjacent layers at the overlapping samples. If the prediction deviation exceeds a preset threshold, a coefficient fine-tuning algorithm is used to correct the regression equation coefficients of adjacent layers, so that the prediction values ​​between layers transition continuously. The local mapping relationships of all coverage layers are concatenated in the order of coverage intervals to form a global nonlinear mapping model: For any input pixel normalized vegetation index value and enhanced vegetation index value, the model first determines its corresponding coverage layer, finds the approximate coverage interval of the pixel in the sample set, and then calls the local mapping relationship of the corresponding layer to calculate the vegetation coverage prediction value. The output result is the vegetation coverage prediction value of the pixel.

[0137] Based on the normalized vegetation index (NDI) sequence, enhanced vegetation index (EMI) sequence, and global nonlinear mapping model at each time point, a pixel-by-pixel analysis algorithm is used to extract the corresponding values ​​of each pixel in the two index sequences within the sub-regions of real cultivated land. These values ​​are then synchronously input into the global nonlinear mapping model to obtain the initial coverage value of the pixel. A spatial consistency check algorithm is used to compare the initial coverage values ​​of adjacent pixels. If the difference exceeds a preset consistency threshold, interpolation correction is performed based on the coverage values ​​of surrounding pixels to ensure the continuity of spatial distribution. All corrected pixel coverage values ​​are arranged according to their spatial coordinates to form a cultivated land vegetation coverage that is spatially heterogeneous and uses pixels as the basic unit.

[0138] Based on farmland vegetation cover, an interval partitioning algorithm is used to divide the cover value into several continuous intervals at preset intervals, where the interval size is determined according to the standard deviation of the sample cover. A weighted statistical algorithm is used to assign weights to pixels within each interval based on their distance from the center of the sub-region; the closer the distance, the higher the weight, reflecting the representativeness of the core area. The proportion of weighted pixels to the total number of weighted pixels is calculated as the interval probability value. A function optimization algorithm is used to fit the cover distribution probability function using the maximum likelihood estimation method, with the cover interval as the horizontal axis and the interval probability value as the vertical axis. This function needs to be verified for goodness of fit through a chi-square test to ensure it accurately reflects the spatial probability distribution characteristics of the cover.

[0139] Based on the absolute value of the change in the vegetation cover of cultivated land and the time characteristics of the same pixel location at continuous time points, all time points are sorted according to the acquisition time through a time axis alignment algorithm to ensure that the interval between adjacent time points is uniform.

[0140] A pixel-level temporal correlation algorithm is used to extract the coverage value of each pixel at consecutive time points, forming a coverage time series for that pixel. The absolute value of the corresponding temporal feature change value is used as the weight of each time point in the series (a preset basic weight is used when the value is zero). A dynamic trend fitting algorithm is used to fit the trend of the coverage time series of each pixel using the weighted least squares method, obtaining the coverage time sub-function of that pixel. A spatial integration algorithm is used to concatenate the coverage time sub-functions of all pixels according to their spatial coordinates, forming a binary function containing spatial location parameters and time parameters, namely the cultivated land vegetation coverage time distribution function. This function can output the coverage prediction value at any time point and any spatial location.

[0141] The first trend unit is used to invert the moisture content distribution state function of the corresponding real farmland sub-region based on the farmland vegetation coverage and the corresponding coverage distribution probability function of each real farmland sub-region, combined with the surface temperature and thermal inertia model. At the same time, based on the moisture content distribution state function of the real farmland sub-region combined with the farmland vegetation coverage time distribution function, the moisture content time fluctuation function of the real farmland sub-region at continuous time points is obtained.

[0142] It should be further explained that one implementation of the moisture content distribution state function of the actual cultivated land sub-regions in this embodiment includes:

[0143] Based on remote sensing images at corresponding time points, thermal radiation values ​​of each pixel within the sub-regions of real cultivated land are extracted using a thermal infrared band analysis algorithm. These thermal radiation values ​​are then converted into surface temperature values, forming surface temperature distribution data that includes the spatial coordinates of all pixels and their corresponding temperature values. Simultaneously, based on the vegetation cover of the cultivated land, vegetation cover correction is applied to the surface temperature distribution data. Pixels with high coverage are assigned a vegetation shading coefficient, where higher coverage results in a larger vegetation shading coefficient, reducing the interference of vegetation on soil temperature and obtaining corrected surface temperature distribution data.

[0144] Based on the corrected surface temperature distribution data and the thermal inertia model, which reflects the relationship between soil thermal inertia and diurnal temperature range and thermal conductivity, the larger the thermal inertia, the slower the soil temperature changes with the environment. The diurnal temperature range is obtained for each pixel in the sub-region through the diurnal temperature difference algorithm, which is the difference between the highest daytime temperature and the lowest nighttime temperature. Combined with the shortwave albedo data of remote sensing imagery, the thermal inertia model is used to calculate the thermal inertia value of each pixel, forming the spatial distribution data of thermal inertia.

[0145] Based on farmland samples with known moisture content, including samples with different vegetation cover and soil types, the thermal inertia value and actual moisture content value of the samples are obtained by synchronous measurement to form a thermal inertia-moisture content sample set. The data in the sample set are fitted by regression analysis algorithm to obtain the basic mapping relationship between thermal inertia value and moisture content value. The thermal inertia increases with the increase of moisture content because the heat capacity of water is higher than that of dry soil.

[0146] Based on the coverage distribution probability function and cultivated land vegetation coverage, a stratified correction algorithm is used to divide the sub-region into several layers according to coverage intervals, including high coverage, medium coverage, and low coverage. Each layer corresponds to a probability interval in the coverage distribution probability function. For each layer sample, the mapping relationship between thermal inertia and water content is refitted. A vegetation transpiration correction coefficient is assigned to the high coverage layer to reduce the sensitivity of thermal inertia to water content, and a soil bareness correction coefficient is assigned to the low coverage layer to increase the sensitivity of thermal inertia to water content, thus obtaining the local mapping relationship for each layer.

[0147] Based on the spatial distribution data of thermal inertia, the local mapping relationship of each layer, and the coverage distribution probability function, a spatial interpolation algorithm is used to determine the layer to which each pixel belongs based on its coverage value. The local mapping relationship of the corresponding layer is then invoked to convert the thermal inertia value into an initial moisture content value. A probability weighting algorithm is then used to adjust the initial moisture content value by combining the probability value corresponding to the coverage of the pixel in the coverage distribution probability function. The higher the probability, the smaller the adjustment range, ensuring that the result conforms to the overall distribution characteristics.

[0148] The adjusted moisture content values ​​of all pixels are integrated according to spatial coordinates to form a moisture content distribution state function that describes the spatial distribution law of moisture content in a sub-region. The input of this function is the spatial coordinates of the pixels, and the output is the corresponding predicted moisture content value.

[0149] It should be further explained that one implementation method of the water content time fluctuation function corresponding to the real cultivated land segmentation sub-regions at continuous time points in this embodiment includes:

[0150] Based on the collection times of continuous time points and the spatial coordinates of the actual cultivated land sub-regions, a time axis calibration algorithm is used to arrange all time points in chronological order to ensure that the intervals between adjacent time points are uniform. If there are differences in intervals, interpolation is performed to complete the data, forming a time series axis with equal time intervals, so that the moisture content data at different time points are comparable in time.

[0151] Based on the spatial coordinates of each pixel within the sub-region of real cultivated land and the state function of moisture content distribution at each time point, a spatiotemporal indexing algorithm is used to extract the moisture content value of the pixel at each time point, and arrange them in order according to the time series axis to form the moisture content time series of the pixel. Since some pixels at certain time points are not included in the sub-region, the missing moisture content values ​​in the time series are filled by interpolation based on the moisture content values ​​of adjacent time points and the predicted value of cultivated land vegetation cover time distribution function of the pixel, so as to ensure the continuity of the series.

[0152] Based on the temporal distribution function of cultivated land vegetation cover, the temporal change rate of cover for each pixel at consecutive time points is extracted to reflect the increase or decrease of cover over time. Through a conversion algorithm, the temporal change rate of cover is converted into a weight value of the water content time series. The larger the temporal change rate of cover, the higher the weight, reflecting the significant impact of vegetation change on water content fluctuation.

[0153] Based on the water content time series of each pixel and its corresponding weight value, a weighted time series analysis algorithm is used to perform trend decomposition on the water content values ​​in the series, separating the long-term trend component and the short-term fluctuation component.

[0154] Using a function fitting algorithm, with time as the independent variable and moisture content as the dependent variable, and combining the periodic characteristics of short-term fluctuation components, such as daily and weekly variations, the weighted least squares method is used to fit and obtain the time fluctuation sub-function of the moisture content of the pixel.

[0155] The moisture content temporal fluctuation sub-functions of all pixels are correlated and integrated according to their spatial coordinates to form the moisture content temporal fluctuation function of real cultivated land segmented sub-regions at continuous time points. This function can output the predicted value of moisture content fluctuation of any pixel in the sub-region at any time point, reflecting the dynamic change law of moisture content over time.

[0156] Based on the vegetation cover and the corresponding cover distribution probability function of each real cultivated land sub-region, combined with the moisture content distribution state function of each real cultivated land sub-region, the soil fertility distribution function of each real cultivated land sub-region is obtained by inversion. At the same time, based on the soil fertility distribution function and the time distribution function of cultivated land vegetation cover of each real cultivated land sub-region at continuous time points, the time change trend of the first soil fertility at continuous time points is obtained by inversion.

[0157] It should be further explained that one specific implementation of the first soil fertility time change trend in this embodiment is as follows:

[0158] First, based on the spatial variation characteristics of the normalized vegetation index (NVI) sequence and the enhanced vegetation index (EVI) sequence, a vegetation growth status extraction algorithm is used to calculate the difference between the NVI and EVI values ​​at adjacent time points for each pixel. This difference is then divided by the time interval to obtain the rate of change. Based on a threshold range determined in advance through historical data statistics, the rate of change is divided into three discrete levels: high, medium, and low. Based on prior knowledge of agricultural ecology, a correspondence rule between these levels and soil fertility is established: a high rate of change indicates vigorous vegetation growth, corresponding to a high soil fertility level; a medium rate of change indicates stable vegetation growth, corresponding to a medium soil fertility level; and a low rate of change indicates sluggish vegetation growth, corresponding to a low soil fertility level.

[0159] Based on the vegetation cover and its distribution probability function, this study uses a stratified coverage algorithm, employing either equal-interval division or natural breakpoint method, to divide the coverage range into three continuous intervals: high coverage, medium coverage, and low coverage. The integral value of the probability density function within each interval is calculated using the coverage distribution probability function as the interval probability value. This probability value is then normalized to a sum of one, yielding the weight coefficients for each interval. For all pixels within each coverage interval, a multiple linear regression analysis algorithm is used, employing vegetation growth status level as the classification independent variable, water content distribution state function value as the continuous independent variable, and historical measured soil fertility value as the dependent variable. This establishes a quantitative regression model between soil fertility value and the two types of independent variables within each coverage interval.

[0160] Subsequently, a spatial fusion algorithm is used to spatially stitch together the soil fertility estimates calculated by the quantitative regression model for each coverage interval according to pixel coordinates to form a preliminary soil fertility distribution map. A Gaussian filtering algorithm or a mean filtering algorithm is used to perform spatial convolution operation on the distribution map. By adjusting the size of the filtering kernel, the smoothness is controlled to eliminate the spatial discontinuity caused by the interval division, and finally a spatially continuous soil fertility distribution function is generated.

[0161] When retrieving the temporal change trend of soil fertility at continuous time points, based on the soil fertility distribution function at continuous time points, a time series extraction algorithm is used to extract the soil fertility values ​​at each time point according to pixel coordinates, and arrange them in chronological order to form a fertility time series for a single pixel. At the same time, based on the temporal distribution function of cultivated land vegetation cover, the difference between the cover value at each time point and the previous time point is calculated and divided by the time interval to obtain the cover time change rate. The range normalization method is used to transform the change rate to the range of 0 to 1 as a weighting factor.

[0162] The weighted trend analysis algorithm is used to multiply the value of each time point in the fertility time series of each pixel by the corresponding weight factor to obtain a new weighted time series. The seasonal decomposition algorithm is used to decompose the new weighted time series into a superposition of long-term trend component, seasonal component and residual component. Based on the fixed step size sliding window algorithm, the difference between the trend value of the next time point and the trend value of the previous time point within the window is calculated with the window size as the unit length to obtain the initial fertility change series.

[0163] Combining the trend characteristics of the temporal distribution function of cultivated land vegetation cover, the slope value of the cover sequence is calculated as a long-term change characteristic through linear regression analysis. A change correction algorithm is adopted, which multiplies the initial change by a decay coefficient less than 1 when the slope is negative, multiplies the initial change by an enhancement coefficient greater than 1 when the slope is positive, and keeps the initial change unchanged when the slope is close to zero. Finally, through a spatial aggregation algorithm, the changes of all pixels after correction are arranged in chronological order, and the average value of the changes of all pixels at each time point is calculated to form the first soil fertility temporal change trend that represents the overall fertility change trend of the region.

[0164] The second trend unit, based on the spatial change values ​​of normalized vegetation index and enhanced vegetation index combined with the crop phenological extraction algorithm, obtains the crop phenological nodes corresponding to the real cultivated land segmented sub-regions. Based on the temporal feature change values ​​of the crop phenological nodes corresponding to the real cultivated land segmented sub-regions, it performs phenological anomaly identification, obtains the corresponding crop phenological anomaly identification feature space, and inverts the second soil fertility temporal change trend of the corresponding real cultivated land segmented sub-region based on the corresponding crop phenological anomaly identification feature space.

[0165] It should be further explained that the process of obtaining the crop phenological nodes corresponding to the actual cultivated land sub-regions in this embodiment includes:

[0166] Based on the spatial variation values ​​of the Normalized Difference Vegetation Index (NDC) and the Enhanced Difference Vegetation Index (EDI), a time series reconstruction algorithm is used to extract the NDC spatial variation values ​​of each pixel at consecutive time points and arrange them in chronological order to form the NDC spatial variation time series for that pixel. Simultaneously, the EDI spatial variation values ​​of each pixel at consecutive time points are extracted and arranged in chronological order to form the EDI spatial variation time series for that pixel. A region aggregation algorithm is used to align the NDC spatial variation time series of all pixels within the sub-regions of the real cultivated land, and the arithmetic mean of the corresponding values ​​of all pixels at each time point is calculated to obtain the sub-region-level NDC spatial variation time series. The same method is used to obtain the sub-region-level EDI spatial variation time series.

[0167] Based on two types of exponential spatial variation time series at the sub-region level, a moving average filtering algorithm is used. A fixed-length window of three time points is employed. For each non-endpoint time point in the sequence, the arithmetic mean of the values ​​from the preceding, current, and following time points is calculated, and this mean is assigned to the current time point as the smoothed value. The original values ​​of the time points at both ends of the sequence are retained. A trend consistency verification algorithm is used to calculate the difference in linear regression slopes of the smoothed and unsmoothed sequences over the overall time range. If the absolute value of the slope difference exceeds a preset threshold, the window length is reduced by one time point, and the moving average calculation is repeated until the slope difference meets the requirements. Finally, the denoised normalized vegetation index spatial variation smoothed sequence and the enhanced vegetation index spatial variation smoothed sequence are obtained.

[0168] Based on the two types of smoothed sequences after denoising, the sequence slope value at each time point is calculated using a numerical differentiation algorithm; the time point where the slope value changes from negative to positive is identified as the rising inflection point, and the time point where the slope value changes from positive to negative is identified as the falling inflection point; the time point with the largest absolute value of the slope is identified as the extreme point; through a multi-sequence association verification algorithm, the positions in the normalized vegetation index spatial change smoothed sequence and the enhanced vegetation index spatial change smoothed sequence that are both identified as feature points of the same type at the same time point are selected and marked as candidate phenological points.

[0169] Based on the candidate phenological point and crop phenological stage extraction algorithm, the earliest rising inflection point is classified as the greening stage node, the first extreme point in the rising stage is classified as the jointing stage node, the time point corresponding to the global maximum value of the sequence is classified as the grain-filling stage node, and the first extreme point in the falling stage is classified as the maturity stage node. The phenological stage time sequence verification algorithm calculates the time interval between the greening stage node and the jointing stage node. If the number of days is not within the preset range specified by the typical growth cycle of the crop, the extreme point with the second largest absolute value of the slope in the rising stage is reselected as the jointing stage node, and the time interval is re-verified until the requirements are met.

[0170] Based on the adjusted phenological nodes, a pixel-level phenological node extraction algorithm is used to obtain the time points of the greening-up stage, jointing stage, grain-filling stage, and maturity stage corresponding to each pixel in the sub-region. A spatiotemporal consistency verification algorithm is used to calculate the absolute value of the time deviation between the phenological time point of each pixel and the corresponding phenological node in the sub-region. Pixels with time deviations within a preset tolerance number of days are retained, and for pixels exceeding the tolerance number of days, their phenological time points are adjusted to the time of the phenological node in the sub-region. Finally, the greening-up stage, jointing stage, grain-filling stage, and maturity stage nodes adjusted for spatiotemporal consistency are integrated in chronological order to form a set of crop phenological nodes corresponding to the sub-regions of real cultivated land.

[0171] It should be further explained that the process of obtaining the temporal change trend of the second soil fertility corresponding to the actual cultivated land sub-regions in this embodiment includes:

[0172] Based on crop phenological nodes corresponding to real cultivated land sub-regions with more than five consecutive historical time points, a time series statistical algorithm is used to calculate the arithmetic mean of the dates of the greening, jointing, grain-filling, and maturity stages within the same agricultural year as the average time value, and the standard deviation of these date values ​​is calculated as the time fluctuation range. Through a multi-feature extraction algorithm, three key feature parameters are extracted from the spatial variation value sequence of the normalized vegetation index of the same historical period: peak value (the maximum value of the sequence), peak duration (the number of consecutive time points where the value remains above 80% of the peak value), and curve rise rate (the average daily change from the growth start point to the peak point). At the same time, the same three feature parameters are extracted from the spatial variation value sequence of the enhanced vegetation index, forming a dual phenological period index feature benchmark set that includes both time features and index features.

[0173] The specific dates of each phenological node are converted into annual day numbers using a date conversion algorithm, i.e., the day of the year. The original time deviation value is obtained by subtracting the current day number from the average time day number of the corresponding phenological nodes in the benchmark set using a difference calculation algorithm. The weighted phenological time deviation value is obtained by converting the time feature change value of each pixel position to a weight coefficient in the range of 0 to 1 using a minimum-maximum normalization method and multiplying the weight coefficient with the original time deviation value.

[0174] The weighted phenological time deviation value is compared with the historical fluctuation range of the corresponding node in the phenological period index feature benchmark set. The upper limit of the historical fluctuation range is the average time value plus twice the standard deviation, and the lower limit is the average time value minus twice the standard deviation. Deviation values ​​exceeding the upper limit are marked as delayed anomalies, and deviation values ​​below the lower limit are marked as premature anomalies. Through a feature comparison algorithm, the relative percentage deviation between the peak value of the current normalized vegetation index spatial change value and the peak value of the corresponding phenological period in the benchmark set is calculated. When the relative deviation percentage is lower than -30%, it is marked as a vitality deficiency anomaly. Through a spatial aggregation algorithm, all anomaly markers are integrated according to phenological period type and pixel coordinates to form a crop phenological period anomaly identification feature space that includes anomaly type, geographic coordinates, and quantified deviation degree.

[0175] Based on the feature space of crop phenological anomalies and historical soil fertility change sample data, a quantitative relationship model between phenological anomalies and the magnitude of soil fertility changes is established using a multiple linear regression algorithm. The number of days of delayed anomaly in the greening stage is taken as continuous independent variable one; the product of the number of days of delayed anomaly in the jointing stage and the degree of insufficient vitality is taken as continuous independent variable two; the weighted sum of the peak anomaly degree in the grain-filling stage and the number of days of delayed anomaly in the maturity stage is taken as continuous independent variable three; and the magnitude of soil fertility change is taken as continuous dependent variable. The regression coefficients are solved using the least squares method. Weight coefficients are assigned to different anomaly types according to the absolute value of the regression coefficients to form a weighted phenological-fertility correlation model.

[0176] The algorithm for calculating changes in soil fertility is used to input the anomaly types and deviation values ​​from the crop phenological anomaly identification feature space into a weighted phenological-fertility correlation model to calculate the soil fertility change for each phenological node. A time series construction algorithm is used to connect the fertility changes of each node in phenological time sequence to form an initial change sequence. A time series smoothing algorithm is then used to process the initial change sequence with a moving average of three time points within a window size to eliminate short-term fluctuations and extract long-term trend components. A trend verification algorithm is used to calculate the slope values ​​of the spatial change trends of the normalized vegetation index and the enhanced vegetation index. When the slope of this slope is opposite in sign to the slope of the fertility change trend, a weighted adjustment is performed on the fertility change trend, with the adjustment weight being the absolute value of the ratio of the two slopes. Finally, a second temporal change trend of soil fertility for the sub-regions of the actual cultivated land is formed.

[0177] The integrated trend unit is used to time-align the first soil fertility time change trend with the second soil fertility time change trend, and then use a weighted average algorithm to combine the first soil fertility time change trend with the second soil fertility time change trend within the aligned time period to obtain the ecological degradation fluctuation value of the real cultivated land segmented sub-region.

[0178] This process, through deep integration of multi-source remote sensing data and crop ecological parameter data, constructs a comprehensive monitoring system encompassing vegetation dynamics, soil moisture, and fertility changes, significantly improving the comprehensiveness and accuracy of farmland ecological status assessment. Its core advantage lies in utilizing the spatial and temporal variation characteristics of normalized and enhanced vegetation indices. This not only achieves high-precision inversion of vegetation cover but also characterizes its spatial heterogeneity through probability distribution functions, laying a reliable foundation for subsequent analysis. Combining surface temperature and thermal inertia models to invert water content distribution fully considers vegetation shading effects and soil thermal characteristics, making water status assessment more closely reflect actual conditions. Furthermore, by fusing the cover time function with water content fluctuations, dynamic tracking of soil fertility is achieved. The first soil fertility temporal change trend reveals fertility changes from the perspective of vegetation growth response and water stress, while the second soil fertility temporal change trend introduces crop physiological time constraints through phenological anomaly identification, effectively capturing abnormal growth rhythms caused by cultivation activities or environmental stress. Finally, the two types of fertility trends were time-aligned and weighted to integrate them, taking into account both the direct coupling effect of vegetation and water and the indirect indicative role of phenological health. This resulted in ecological degradation fluctuation values ​​that have both physical and ecological significance, and can sensitively reflect the gradual degradation or sudden deterioration of arable land ecosystems, providing a quantitative basis for sustainable management of regional arable land resources and ecological restoration decisions.

[0179] The early warning module provides real-time warnings based on fluctuations in farmland area or ecological degradation in actual farmland sub-regions, combined with corresponding preset warning values. It also automatically generates warning information and visualizes the spatial location, type, and degree of change in the actual farmland sub-regions through a GIS platform. It should be further noted that the module in this embodiment includes an early warning discrimination unit and a real-time annotation unit.

[0180] The early warning judgment unit, based on the fluctuation values ​​of cultivated land area and ecological degradation at consecutive time points, determines that there is an anomaly in the sub-region of real cultivated land division at the corresponding consecutive adjacent time points if either the fluctuation value of cultivated land area or the fluctuation value of ecological degradation at two adjacent time points does not meet the corresponding early warning value.

[0181] The real-time annotation unit determines the spatial location and area of ​​the changed region based on the absolute value of the temporal characteristic change of adjacent real cultivated land sub-regions with anomalies at the same pixel location point under continuous time points, the fluctuation value of the first cultivated land area within the corresponding time period, and the corresponding pixel location point. Simultaneously, based on the normalized difference vegetation index (NDVI) and enhanced vegetation index (EGI) within the changed region and their corresponding spatial change values, it determines the type of change. Based on the determined spatial location and area of ​​the changed region, combined with the magnitude of the absolute value of the temporal characteristic change within the corresponding continuous time period, an evaluation algorithm is used to obtain the degree of change of cultivated land type within the current changed region. The spatial location, area of ​​change, type of change, and degree of change are then dynamically annotated and displayed in real-time on the real cultivated land sub-regions corresponding to adjacent time points in the spatiotemporal coordinate system using an automatic annotation algorithm.

[0182] It should be further explained that the anomaly detection process in this embodiment includes:

[0183] Based on the fluctuation values ​​of cultivated land area and ecological degradation in sub-regions of actual cultivated land at continuous historical time points, a statistical analysis algorithm is used to calculate the historical maximum value and the upper limit of the 95% confidence interval for cultivated land area fluctuation. This upper limit is determined as the warning value for cultivated land area fluctuation; exceeding this value indicates that the area change is outside the normal range. Similarly, the historical maximum value and the upper limit of the 95% confidence interval for ecological degradation fluctuation are calculated, and this upper limit is determined as the warning value for ecological degradation fluctuation. Exceeding this value indicates that the degree of ecological degradation is abnormal.

[0184] When extracting fluctuation values ​​at adjacent time points, based on the sequence of fluctuation values ​​of cultivated land area and the sequence of fluctuation values ​​of ecological degradation at continuous time points, the sequence truncation algorithm is used to select the fluctuation values ​​of cultivated land area corresponding to any two adjacent time points, which are recorded as the current area fluctuation value and the previous area fluctuation value, and the corresponding ecological degradation fluctuation value, which are recorded as the current degradation fluctuation value and the previous degradation fluctuation value, to form a pair of fluctuation values ​​at adjacent time points.

[0185] When performing fluctuation value comparison, the current area fluctuation value is compared with the cultivated land area fluctuation warning value through a threshold comparison algorithm to determine whether the current area fluctuation value is less than or equal to the cultivated land area fluctuation warning value, i.e., the warning value is met. At the same time, the previous area fluctuation value is compared with the cultivated land area fluctuation warning value to determine whether it meets the warning value. Similarly, the current degradation fluctuation value and the previous degradation fluctuation value are compared with the ecological degradation fluctuation warning value to determine whether the warning value is met.

[0186] When an anomaly is determined in a region, a comprehensive judgment is made based on the above comparison results using a logical judgment algorithm: if, in two adjacent time points, the current area fluctuation value does not meet the warning value for cultivated land area fluctuation (i.e., it is greater than the warning value), or the previous area fluctuation value does not meet the warning value for cultivated land area fluctuation, or the current degradation fluctuation value does not meet the warning value for ecological degradation fluctuation, or the previous degradation fluctuation value does not meet the warning value for ecological degradation fluctuation, then it is determined that the actual cultivated land sub-region under the consecutive adjacent time points is abnormal; only when all comparison results meet the corresponding warning values ​​is it determined to be normal.

[0187] It should be further explained that the implementation process of the automatic annotation algorithm in this embodiment includes:

[0188] Based on the identification of anomalies in adjacent real farmland sub-regions at continuous time points, a pixel-level change detection algorithm is used to extract the absolute value of the temporal feature change value of each pixel position at adjacent time points. This absolute value is compared with a preset fluctuation threshold, and all pixels exceeding the threshold are filtered out. A spatial clustering analysis algorithm is used, employing the density-based DBSCAN clustering method. A neighborhood radius parameter and a minimum pixel number parameter are set to cluster pixels that are within the neighborhood radius and meet the minimum pixel number requirement, forming a continuous set of spatial coordinates of the change area. At the same time, isolated pixels that do not meet the density requirements are removed.

[0189] Based on the set of spatial coordinates of the changed region, the actual geographic area corresponding to each pixel is calculated by using a pixel area conversion algorithm according to the spatial resolution parameters of the remote sensing image, i.e., the actual ground size represented by a single pixel.

[0190] The initial changed area is obtained by summing the actual geographical areas of all pixels within the changed area using a regional area statistics algorithm. The area calibration algorithm calculates the area calibration coefficient using the fluctuation value of the first cultivated land area at adjacent time points. Specifically, the fluctuation value of the first cultivated land area is normalized and then multiplied by one as a multiplier, which is then multiplied by the initial changed area to obtain the calibrated actual changed area.

[0191] When determining the type of change, the difference between index values ​​at adjacent time points is calculated based on the normalized vegetation index sequence and enhanced vegetation index sequence of each pixel in the change area through a trend analysis algorithm.

[0192] According to the type determination rule: if the difference between the normalized vegetation index and the difference between the enhanced vegetation index are both less than zero and the spatial variation value is negative, it is determined to be a type of farmland degradation.

[0193] If the difference between the two indices is greater than zero and the spatial change value is positive, it is determined to be a farmland restoration type; if the signs of the difference between the two indices are inconsistent or the absolute value of the spatial change value exceeds the disorder threshold, it is determined to be a farmland use transformation type.

[0194] When assessing the degree of change, a tiered assessment algorithm is used to divide the actual area of ​​change into three levels (small, medium, and large) based on the tertiles, and the absolute value of the time characteristic change value is divided into three levels (low, medium, and high) based on the tertiles. A weighted scoring algorithm is used to sum the actual area of ​​change level with a weight of 60% and the absolute value of the time characteristic change value with a weight of 40% to obtain a comprehensive score. A severity classification algorithm is used to divide the comprehensive score into three severity levels (mild, moderate, and severe) based on the tertile range.

[0195] When performing real-time dynamic annotation, a coordinate transformation algorithm maps the spatial coordinates of the changed area to geographic coordinates in the spatiotemporal coordinate system; a visualization annotation algorithm uses red semi-transparent polygons to annotate areas of farmland degradation, green semi-transparent polygons to annotate areas of farmland restoration, and yellow semi-transparent polygons to annotate areas of farmland utilization transformation; a size mapping algorithm adjusts the display size of the annotated polygons proportionally according to the actual changed area value; and a color depth mapping algorithm adjusts the color transparency of the annotated polygons according to the severity level, with the highest transparency for severe levels and the lowest transparency for mild levels, forming a dynamically updated visualization annotation display.

[0196] This process, through the construction of multi-dimensional fluctuation indicators and intelligent discrimination mechanisms, achieves efficient and accurate early warning of abnormal changes in arable land resources, significantly improving the real-time performance and decision support capabilities of arable land dynamic monitoring. Its core value lies in the comprehensive utilization of both area fluctuation and ecological degradation fluctuation indicators, capturing both rapid changes in arable land quantity and revealing the gradual degradation of ecosystem quality, overcoming the one-sidedness of single-indicator early warning and making anomaly judgment more comprehensive and reliable. Early warning thresholds are determined through historical data statistics, fully reflecting regional specificity and the normal fluctuation range, avoiding false alarms and missed alarms. Furthermore, the automatic labeling process deeply integrates pixel-level change detection, spatial clustering, and multi-source feature analysis, accurately locating the spatial extent of changed areas and quantifying the changed area. Simultaneously, it intelligently discriminates change types based on vegetation index change trends and spatial characteristics, clearly distinguishing between different scenarios such as degradation, restoration, and utilization transformation. Severity assessment comprehensively considers both the changed area and change intensity, employing a weighted scoring mechanism to ensure the scientific rigor and intuitiveness of the grading results. Finally, dynamic visualization annotations are performed through a GIS platform, transforming abstract data fluctuations into concrete spatial graphics. Color, transparency, and size are used to map the type and severity of changes, providing managers with an extremely intuitive and information-rich monitoring view.

[0197] Example 2:

[0198] Please see Figure 2 Another embodiment of the present invention provides a method for monitoring land ecological conditions based on remote sensing data, comprising the following steps:

[0199] S1. Obtain remote sensing image sequences of the target area and analyze them to obtain initial farmland segmentation sub-region sequences and edge segmentation feature sets;

[0200] S2. Based on the initial cultivated land segmentation sub-region sequence and edge segmentation feature set, obtain the remote sensing image feature space;

[0201] S3. Based on the spatial features of remote sensing images and the initial farmland sub-regions, determine the area of ​​the actual farmland sub-regions, and based on the determined changes in the area of ​​the actual farmland sub-regions, determine the farmland area fluctuation value.

[0202] S4. Based on the spatial characteristics of remote sensing images, combined with surface temperature and inversion, the vegetation cover of cultivated land is obtained, and the ecological degradation fluctuation value of the real cultivated land segmented sub-regions is obtained.

[0203] S5. Real-time early warning is issued based on the fluctuation value of cultivated land area or the fluctuation value of ecological degradation in the actual cultivated land sub-regions, combined with the corresponding preset early warning value. At the same time, early warning information is automatically generated, and the spatial location, type and degree of change of the actual cultivated land sub-regions are visualized and marked through the GIS platform.

[0204] The embodiments of the present invention have been described above with reference to the accompanying drawings. However, the present invention is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments under the guidance of the present invention without departing from the spirit and scope of the claims. All of these variations are within the protection scope of the present invention.

Claims

1. A land ecological status monitoring system based on remote sensing data, characterized in that, include: The module includes a preprocessing module, a feature extraction module, a region determination module, an inversion module, and an early warning module. The preprocessing module is used to acquire remote sensing image sequences of the target area and analyze and acquire initial farmland segmentation sub-region sequences and edge segmentation feature sets; The feature extraction module obtains a remote sensing image feature space based on the initial cultivated land segmentation sub-region sequence and edge segmentation feature set. The remote sensing image feature space includes spatial feature sequences and spatial feature fluctuation values ​​for any N edge extension directions of the remote sensing image. The spatial feature sequences for any N edge extension directions include normalized vegetation index sequences and enhanced vegetation index sequences for the corresponding extension directions. The spatial feature fluctuation values ​​include the spatial variation values ​​of the normalized vegetation index and enhanced vegetation index for the remote sensing image at each time point. The region determination module determines the area of ​​the actual cultivated land sub-region based on the spatial features of remote sensing images combined with the initial cultivated land sub-regions, and determines the cultivated land area fluctuation value based on the determined change value of the actual cultivated land sub-region area. The inversion module, based on the spatial features of remote sensing images, combines surface temperature with inversion to obtain the vegetation cover of cultivated land, and obtains the ecological degradation fluctuation value of the real cultivated land sub-regions. The early warning module provides real-time early warnings based on the fluctuation values ​​of cultivated land area or ecological degradation in the actual cultivated land sub-regions, combined with corresponding preset early warning values. It also automatically generates early warning information and visualizes and marks the spatial location, type, and degree of change of the actual cultivated land sub-regions through a GIS platform.

2. The land ecological status monitoring system based on remote sensing data as described in claim 1, characterized in that, The preprocessing module includes an acquisition unit and an image segmentation unit; The acquisition unit is used to acquire a sequence of remote sensing images of the target area over a continuous preset time period, and to perform radiometric and geometric correction preprocessing on the remote sensing image sequence to obtain a corrected remote sensing image sequence. The image segmentation unit is used to perform initial region identification based on the corrected remote sensing image sequence and the pre-trained YOLOv5 model, to obtain the initial region identification range and corresponding region type of the remote sensing images collected at different time points, and to obtain the edge segmentation feature set corresponding to the initial cultivated land segmentation sub-region at each time point based on the initial region identification range and the image segmentation algorithm; the region type includes cultivated land, construction land, forest land, pond, and abandoned land.

3. The land ecological status monitoring system based on remote sensing data as described in claim 2, characterized in that, The feature extraction module includes a center determination unit and a feature extraction unit; The center determination unit obtains a first center point based on the pixel position points in the edge segmentation feature set corresponding to the initial cultivated land segmentation sub-region at each time point, combined with the minimum bounding rectangle, and obtains the centroid position as a second center point based on the edge pixel position points; and determines the true center position point of the initial cultivated land segmentation sub-region at each time point based on the average of the first center point and the second center point. Two-dimensional coordinates are constructed based on the initial farmland sub-regions of the real center location at each time point, and the line connecting all aligned real center locations at consecutive time points is used as a three-dimensional time axis to construct a spatiotemporal coordinate system. The initial farmland sub-regions of the remote sensing images collected at each time point are embedded into the time point corresponding to the spatiotemporal coordinate system. The feature extraction unit, based on the initial cultivated land sub-region at each time point, takes the real center point of the corresponding initial cultivated land sub-region as the starting point and performs feature extraction along any N edge extension directions in the initial cultivated land sub-region at the corresponding time point to obtain the spatial feature sequence and spatial feature fluctuation value of any N edge extension directions at the current time point.

4. The land ecological status monitoring system based on remote sensing data as described in claim 3, characterized in that, The region determination module includes a region determination unit and a region area unit; The region determination unit, based on the spatial feature sequence and spatial feature fluctuation value of any N edge extension directions at the current time point, combined with the region type fluctuation threshold, if the spatial feature fluctuation value of M consecutive pixel positions in any edge extension direction does not meet the region type fluctuation threshold, then the pixel position point corresponding to the largest spatial feature fluctuation value is taken as the edge point of the current edge extension direction, and this is repeated N times to determine the second edge position point set. The area unit is obtained by combining the second edge position point set of the initial cultivated land sub-region at each time point with the corresponding initial cultivated land sub-region using a region integration algorithm to obtain the area of ​​the actual cultivated land sub-region at each time point. At the same time, taking the pixel point corresponding to the largest actual cultivated land sub-region as the starting point, along the positive and negative axes of the three-dimensional time axis, the absolute value of the time feature change value of the area of ​​adjacent actual cultivated land sub-regions at the same pixel position point at consecutive time points is extracted.

5. The land ecological status monitoring system based on remote sensing data as described in claim 4, characterized in that, The region determination module further includes a region fluctuation determination unit; the region fluctuation determination unit is used to obtain a first cultivated land area fluctuation value sequence based on the area of ​​the actual cultivated land sub-regions corresponding to adjacent time points, and simultaneously obtain the proportion of the absolute values ​​of the time feature changes of the same pixel location points that do not meet the region type fluctuation threshold based on the absolute value of the time feature changes combined with the region type fluctuation threshold, as a second cultivated land area fluctuation value sequence; and obtain the cultivated land area fluctuation value at adjacent time points based on the average of the first cultivated land area fluctuation value and the second cultivated land area fluctuation value at corresponding adjacent time points in the first cultivated land area fluctuation value sequence and the second cultivated land area fluctuation value sequence.

6. The land ecological status monitoring system based on remote sensing data as described in claim 5, characterized in that, The inversion module includes a coverage inversion unit; The coverage inversion unit, based on the normalized vegetation index sequence and enhanced vegetation index sequence corresponding to each time point, and the spatial variation values ​​of the normalized vegetation index and enhanced vegetation index at the corresponding time point, inverts to obtain the cultivated land vegetation coverage and the corresponding coverage distribution probability function for each real cultivated land sub-region at the corresponding time point. At the same time, based on the cultivated land vegetation coverage corresponding to each real cultivated land sub-region and the absolute value of the temporal feature change value of the area of ​​adjacent real cultivated land sub-regions at the same pixel position at consecutive time points, it obtains the cultivated land vegetation coverage time distribution function.

7. The land ecological status monitoring system based on remote sensing data as described in claim 6, characterized in that, The inversion module further includes a first trend unit; The first trend unit is used to invert the moisture content distribution state function of the corresponding real farmland sub-region based on the farmland vegetation coverage and the corresponding coverage distribution probability function of each real farmland sub-region, combined with the surface temperature and thermal inertia model. At the same time, based on the moisture content distribution state function of the real farmland sub-region combined with the farmland vegetation coverage time distribution function, the moisture content time fluctuation function of the real farmland sub-region at continuous time points is obtained. Based on the vegetation cover and the corresponding cover distribution probability function of each real cultivated land sub-region, and combined with the moisture content distribution state function of each real cultivated land sub-region, the soil fertility distribution function of each real cultivated land sub-region is obtained by inversion. At the same time, based on the soil fertility distribution function and the cultivated land vegetation cover time distribution function of each real cultivated land sub-region at continuous time points, the first soil fertility time change trend at continuous time points is obtained by inversion.

8. The land ecological status monitoring system based on remote sensing data as described in claim 7, characterized in that, The inversion module also includes a second trend unit and a comprehensive trend unit; The second trend unit, based on the spatial change values ​​of normalized vegetation index and enhanced vegetation index combined with the crop phenological extraction algorithm, obtains the crop phenological nodes corresponding to the real cultivated land segmented sub-regions. Based on the temporal feature change values ​​corresponding to the crop phenological nodes corresponding to the real cultivated land segmented sub-regions, it performs phenological anomaly identification, obtains the corresponding crop phenological anomaly identification feature space, and inverts the second soil fertility temporal change trend of the corresponding real cultivated land segmented sub-region based on the corresponding crop phenological anomaly identification feature space. The integrated trend unit is used to perform time alignment processing on the first soil fertility time change trend and the second soil fertility time change trend, and to obtain the ecological degradation fluctuation value of the real cultivated land segmented sub-region by combining the first soil fertility time change trend and the second soil fertility time change trend within the aligned time period with a weighted average algorithm.

9. The land ecological status monitoring system based on remote sensing data as described in claim 8, characterized in that, The early warning module includes an early warning discrimination unit and a real-time annotation unit; The early warning judgment unit determines that there is an anomaly in the real farmland sub-region at the corresponding consecutive adjacent time points if either the farmland area fluctuation value or the ecological degradation fluctuation value at two adjacent time points does not meet the corresponding early warning value. The real-time annotation unit determines the spatial location and area of ​​the changed region based on the absolute value of the temporal characteristic change of adjacent real cultivated land sub-regions with anomalies at the same pixel location point under continuous time points, the fluctuation value of the first cultivated land area within the corresponding time period, and the corresponding pixel location point. Simultaneously, based on the normalized difference vegetation index and enhanced vegetation index within the changed region and their corresponding spatial change values, it determines the type of change. Based on the determined spatial location and area of ​​the changed region, combined with the magnitude of the absolute value of the temporal characteristic change within the corresponding continuous time period, it obtains the degree of change of cultivated land type within the current changed region through an evaluation algorithm. The spatial location, area of ​​change, type of change, and degree of change are then dynamically annotated and displayed in real-time on the real cultivated land sub-regions corresponding to adjacent time points in the spatiotemporal coordinate system using an automatic annotation algorithm.

10. A method for monitoring land ecological status based on remote sensing data, implemented using the land ecological status monitoring system based on remote sensing data as described in any one of claims 1-9, characterized in that, include: Acquire remote sensing image sequences of the target area and analyze them to obtain initial farmland segmentation sub-region sequences and edge segmentation feature sets; Based on the initial farmland segmentation sub-region sequence and edge segmentation feature set, the remote sensing image feature space is obtained; Based on the spatial features of remote sensing images combined with the initial farmland sub-regions, the area of ​​the actual farmland sub-regions is determined, and the farmland area fluctuation value is determined based on the determined change value of the actual farmland sub-region area. Based on the spatial characteristics of remote sensing images, combined with surface temperature and inversion, the vegetation cover of cultivated land is obtained, and the ecological degradation fluctuation value of the real cultivated land sub-regions is obtained. Real-time early warnings are issued based on the fluctuation values ​​of cultivated land area or ecological degradation in real cultivated land sub-regions, combined with corresponding preset early warning values. At the same time, early warning information is automatically generated, and the spatial location, type, and degree of change of the real cultivated land sub-regions are visualized and marked through the GIS platform.

Images