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

Abstract

The application 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, which comprises: a preprocessing module that acquires a remote sensing image sequence and obtains an initial cultivated land segmentation sub-region and edge features; a feature extraction module that extracts spatial features such as normalized vegetation index and enhanced vegetation index and wave values based on the same, and constructs a space-time coordinate system; a region determination module that determines a true cultivated land area and a cultivated land area wave value; an inversion module that inverts vegetation coverage, water content and soil fertility, and obtains an ecological degradation wave value in combination with abnormal identification of a phenological period; and an early warning module that performs real-time early warning based on the cultivated land area wave value and the ecological degradation wave value, and visually marks the position, type and severity of a changed region through a GIS platform; and the application realizes accurate monitoring and dynamic early warning of land ecological conditions.

Description

TECHNICAL FIELD

[0001] The application 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. BACKGROUND

[0002] As an important ecological system, the service value and ecological condition evaluation of cultivated land is the key to implement cultivated land protection and natural resource management; in the prior art, a cultivated land resource ecosystem service value accounting method based on multi-source remote sensing technology is disclosed in Chinese Patent No. CN118731935B, which extracts crop categories and planting areas by integrating time-series optical images, radar images and field sample data, constructs a yield correction model coupled with terrain factors and a regional adaptive model of equivalent factors, and realizes cultivated land ecosystem service value accounting; a cultivated land ecosystem service evaluation method based on remote sensing monitoring is disclosed in Chinese Patent Application No. CN119761854A, which obtains multi-dimensional cultivated land data based on remote sensing technology, builds an evaluation model after dimension reduction processing, and realizes cultivated land ecological service evaluation through dynamic weight distribution and service coordination relationship optimization. Both of them try to improve the scientificity of cultivated land ecological evaluation through remote sensing technology, but there are still technical limitations in practical application; among them, the cultivated land resource ecosystem service value accounting method based on multi-source remote sensing technology focuses on value accounting, but it does not include real-time dynamic monitoring mechanism, and the integration of ecological parameters with strong timeliness such as soil moisture and crop phenology is insufficient, which is difficult to reflect the instantaneous change of cultivated land ecology; the cultivated land ecosystem service evaluation method based on remote sensing monitoring focuses on ecological service evaluation, relies on a single remote sensing data source, does not fully integrate multi-source heterogeneous data such as radar and optical images, and the evaluation index is not linked with key ecological parameters such as soil fertility and land surface temperature, resulting in incomplete evaluation dimension. SUMMARY

[0003] In view of the deficiencies of the prior art, the present application provides a land ecological condition monitoring system and method based on remote sensing data, which comprises: a preprocessing module that obtains a sequence of remote sensing images and obtains an initial cultivated land segmentation sub-region and edge features; a feature extraction module that extracts spatial features such as normalized vegetation index and enhanced vegetation index and wave values based on the same, and constructs a spatio-temporal coordinate system; a region determination module that determines the true cultivated land area and the cultivated land area fluctuation value; an inversion module that inverts the vegetation coverage, water content and soil fertility, and obtains the ecological degradation fluctuation value in combination with the abnormal identification of the phenological period; and an early warning module that performs real-time early warning based on the cultivated land area fluctuation value and the ecological degradation fluctuation value, and visualizes the location, type and severity of the changed area through a GIS platform; the present application realizes accurate monitoring and dynamic early warning of the land ecological condition.

[0004] To achieve the above-mentioned purpose, the present application provides the following technical solutions:

[0005] The land ecological condition monitoring system based on remote sensing data comprises a preprocessing module, a feature extraction module, a region determination module, an inversion module and a warning module.

[0006] The preprocessing module is used for acquiring a remote sensing image sequence of a target region and analyzing and acquiring an initial cultivated land segmented sub-region sequence and an edge segmentation feature set.

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

[0008] The region determination module determines a real cultivated land segmented sub-region area based on the remote sensing image feature space in combination with the initial cultivated land segmented sub-region, and determines a cultivated land area fluctuation value based on a determined real cultivated land segmented sub-region area change value.

[0009] The inversion module obtains an ecological degradation fluctuation value of the real cultivated land segmented sub-region based on the remote sensing image feature space in combination with a ground surface temperature and an inversion obtained cultivated land vegetation coverage.

[0010] The warning module performs real-time warning based on the cultivated land area fluctuation value or the ecological degradation fluctuation value of the real cultivated land segmented sub-region in combination with a corresponding preset warning value, simultaneously generates warning information, and visualizes and labels a change region spatial position, a change type and a change degree of the real cultivated land segmented sub-region through a GIS platform.

[0011] Specifically, the preprocessing module comprises a collection unit and an image segmentation unit.

[0012] The collection unit is used for acquiring a remote sensing image sequence of a target region under a continuous preset time length, performing radiation correction and geometric correction preprocessing on the remote sensing image sequence, and obtaining a corrected remote sensing image sequence.

[0013] The image segmentation unit is used for performing region initial identification on the corrected remote sensing image sequence in combination with a pre-trained YOLOv5 model, acquiring an initial region identification range and a corresponding region type of the remote sensing image collected at different time points, and obtaining an edge segmentation feature set corresponding to the initial cultivated land segmented sub-region at each time point based on the initial region identification range in combination with an image segmentation algorithm; the region type comprises cultivated land, construction land, woodland, pond and abandoned land.

[0014] Specifically, the feature extraction module comprises a center determination unit and a feature extraction unit.

[0015] a center determination unit, which obtains a first center point based on a pixel position point in an edge segmentation feature set corresponding to the initial farmland segmentation sub-region at each time point and a minimum circumscribed rectangle, and obtains a second center point as a centroid position based on the edge pixel position point; and determines a real center position point of the initial farmland segmentation sub-region at each time point based on a mean value of the first center point and the second center point;

[0016] a two-dimensional coordinate is constructed based on the real center position point of each time point in the initial farmland segmentation sub-region at the corresponding time, and a three-dimensional time axis is constructed by connecting all the aligned real center position points at continuous time points, to obtain a space-time coordinate system, and the initial farmland segmentation sub-region of the remote sensing image collected at each time point is embedded into the space-time coordinate system at the corresponding time point;

[0017] a feature extraction unit, which extracts features along any N edge extension directions of the initial farmland segmentation sub-region at the corresponding time point based on the real center position point of the initial farmland segmentation sub-region at each time point, to obtain a spatial feature sequence and a 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 the spatial feature fluctuation value of any N edge extension directions at the current time point, in combination with a region type fluctuation threshold value, when the spatial feature fluctuation value of any edge extension direction at continuous M pixel position points does not satisfy the region type fluctuation threshold value, takes the pixel position point corresponding to the maximum spatial feature fluctuation value as an edge point of the current edge extension direction, and repeats N times to determine a second edge position point set.

[0020] The region area unit, based on the second edge position point set of the initial farmland segmentation sub-region at each time point and the corresponding initial farmland segmentation sub-region in combination with a region integral algorithm, obtains the area of the real farmland segmentation sub-region at each time point, and simultaneously takes each position pixel point of the real farmland segmentation sub-region with the largest area as a starting point to extract the absolute value of the time feature change value of the adjacent real farmland segmentation sub-region area at the same pixel position point at continuous time points along the positive and negative axes of the three-dimensional time axis.

[0021] Specifically, the area determination module comprises an area fluctuation determination unit; the area fluctuation determination unit is configured to: obtain a first cultivated land area fluctuation value sequence according to the areas of the real cultivated land sub-regions corresponding to adjacent time points; obtain a proportion of absolute values of the time feature change values of the same pixel position points that do not satisfy the area type fluctuation threshold, as a second cultivated land area fluctuation value sequence, based on the absolute values of the time feature change values in combination with the area type fluctuation threshold; and obtain the cultivated land area fluctuation values of the adjacent time points based on the mean values of the first cultivated land area fluctuation values and the second cultivated land area fluctuation values of the corresponding adjacent time points in the first cultivated land area fluctuation value sequence and the second cultivated land area fluctuation value sequence.

[0022] Specifically, the inversion module comprises a coverage inversion unit.

[0023] The coverage inversion unit is configured to: obtain the cultivated land vegetation coverage corresponding to each real cultivated land sub-region and the coverage distribution probability function of the corresponding time point based on the normalized vegetation index sequence and the enhanced vegetation index sequence corresponding to each time point and the normalized vegetation index spatial change value and the enhanced vegetation index spatial change value corresponding to the corresponding time point; and obtain the cultivated land vegetation coverage time distribution function based on the cultivated land vegetation coverage corresponding to each real cultivated land sub-region and the absolute values of the time feature change values of the adjacent real cultivated land sub-region areas at the same pixel position point at the continuous time points.

[0024] Specifically, the inversion module further comprises a first change trend unit.

[0025] The first change trend unit is configured to: obtain the moisture content distribution state function of the corresponding real cultivated land sub-region based on the cultivated land vegetation coverage corresponding to each real cultivated land sub-region and the coverage distribution probability function in combination with the land surface temperature and the thermal inertia model; and obtain the moisture content time fluctuation function of the real cultivated land sub-region at the continuous time points based on the moisture content distribution state function of the real cultivated land sub-region and the cultivated land vegetation coverage time distribution function.

[0026] The inversion module is configured to: obtain the soil fertility distribution function of each real cultivated land sub-region based on the cultivated land vegetation coverage corresponding to each real cultivated land sub-region and the coverage distribution probability function in combination with the moisture content distribution state function of each real cultivated land sub-region; and obtain the first soil fertility time change trend at the continuous time points based on the soil fertility distribution function of each real cultivated land sub-region at the continuous time points and the cultivated land vegetation coverage time distribution function.

[0027] Specifically, the inversion module further comprises a second change trend unit and a comprehensive trend unit.

[0028] The second change trend unit obtains a crop phenological node corresponding to the real farmland segmentation sub-region based on the normalized vegetation index spatial change value and the enhanced vegetation index spatial change value combined with the crop phenological period extraction algorithm, performs phenological period anomaly recognition based on a time characteristic change value corresponding to the crop phenological node corresponding to the real farmland segmentation sub-region, obtains a corresponding crop phenological period anomaly recognition feature space, and inversely obtains a second soil fertility time change trend of the corresponding real farmland segmentation sub-region based on the corresponding crop phenological period anomaly recognition feature space.

[0029] The comprehensive trend unit is configured to perform time alignment processing on the first soil fertility time change trend and the second soil fertility time change trend, and obtain an ecological degradation fluctuation value of the real farmland segmentation sub-region by using a combined weighted average algorithm of the first soil fertility time change trend and the second soil fertility time change trend in the aligned time period.

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

[0031] The early warning discrimination unit is configured to, based on the farmland area fluctuation value and the ecological degradation fluctuation value at consecutive time points, determine that there is an anomaly in the real farmland segmentation sub-region at the corresponding consecutive adjacent time points if one of the corresponding farmland area fluctuation values or the ecological degradation fluctuation values at adjacent two time points does not satisfy the corresponding early warning value.

[0032] The real-time labeling unit is configured to determine the spatial position and the changed area of the change region based on the absolute value of the time characteristic change value of the area of the adjacent real farmland segmentation sub-region determined to have an anomaly at the consecutive time points at the same pixel position point, the first farmland area fluctuation value in the corresponding time period, and the corresponding pixel position point, determine the type of the change based on the corresponding normalized vegetation index and enhanced vegetation index and the corresponding spatial change value in the change region, obtain the change degree of the farmland type in the current change region by an evaluation algorithm based on the spatial position and the changed area of the change region, the size of the absolute value of the time characteristic change value in the corresponding consecutive time period, and display the spatial position and the changed area, the type of the change, and the change degree in real time and dynamically on the real farmland segmentation sub-region at the corresponding adjacent time points in the spatiotemporal coordinate system by an automatic labeling algorithm.

[0033] The land ecological condition monitoring method based on remote sensing data includes:

[0034] Obtaining a remote sensing image sequence of a target region and analyzing and obtaining an initial farmland segmentation sub-region sequence and an edge segmentation feature set;

[0035] Obtaining a remote sensing image feature space based on the initial farmland segmentation sub-region sequence and the edge segmentation feature set;

[0036] Determine the real cultivated land segmentation sub-region area based on the remote sensing image feature space combined with the initial cultivated land segmentation sub-region, determine the cultivated land area fluctuation value based on the determined real cultivated land segmentation sub-region area change value;

[0037] Based on the remote sensing image feature space combined with the ground surface temperature and the obtained cultivated land vegetation coverage, the ecological degradation fluctuation value of the real cultivated land segmentation sub-region is obtained;

[0038] Based on the real cultivated land segmentation sub-region area fluctuation value or the ecological degradation fluctuation value combined with the corresponding preset warning value, real-time warning is carried out, and warning information is automatically generated, and the change area spatial position, change type and change degree of the real cultivated land segmentation sub-region are visually marked through the GIS platform.

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

[0040] The present application is aimed at the deficiencies of the prior art, and through the construction of a multi-module collaborative remote sensing monitoring system, high-precision dynamic perception and intelligent warning of the ecological condition of cultivated land are realized. The preprocessing module improves the accuracy of initial segmentation of cultivated land through radiation and geometric correction combined with YOLOv5 region recognition. The feature extraction module realizes spatial feature fluctuation quantification based on spatiotemporal coordinate system and multi-direction vegetation index sequence analysis. The area determination module accurately captures the area change of cultivated land through edge point dynamic calibration and area integral algorithm. The inversion module realizes multi-dimensional inversion of soil fertility and ecological degradation trend by fusing vegetation coverage, ground surface temperature, water content and phenological characteristics. The warning module realizes real-time dynamic display of the position, type and severity of the change area by combining GIS visualization and automatic labeling technology. The present application effectively improves the automation, precision and real-time performance of land ecological monitoring, and provides reliable technical support for cultivated land protection and ecological management. BRIEF DESCRIPTION OF DRAWINGS

[0041] Figure 1 Figure 1 is a module diagram of the land ecological condition monitoring system based on remote sensing data according to Embodiment 1 of the present application;

[0042] Figure 2 Figure 2 is a flowchart of the land ecological condition monitoring method based on remote sensing data according to Embodiment 2 of the present application. DETAILED DESCRIPTION

[0043] Embodiment 1:

[0044] Please refer to Figure 1 An embodiment provided by the present application is a land ecological condition monitoring system based on remote sensing data, which comprises a preprocessing module, a feature extraction module, an area determination module, an inversion module and a warning module.

[0045] The preprocessing module is configured to acquire a remote sensing image sequence of a target region and analyze the acquired initial cultivated land segmentation sub-region sequence and edge segmentation feature set. It should be further explained that the preprocessing module in the embodiment includes a collection unit and an image segmentation unit.

[0046] The collection unit is configured to acquire a remote sensing image sequence of the target region within a continuous preset time length, perform radiation correction and geometric correction preprocessing on the remote sensing image sequence, and obtain a corrected remote sensing image sequence.

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

[0048] Based on a digital quantization value of the remote sensing image, the digital quantization value is an integer with a preset number of bits, and the value is within a preset minimum value to a preset maximum value range. Through a sensor calibration algorithm, a gain coefficient provided by a satellite sensor is used, wherein the gain coefficient adopts a preset unit and the value is within a preset gain coefficient range. An offset coefficient is used, which adopts the same preset unit as the gain coefficient and the value is within a preset offset coefficient range. The digital quantization value is converted into a radiation brightness value received by the sensor, and the radiation brightness value is within a preset radiation brightness range, and the accuracy meets a preset radiation brightness accuracy standard.

[0049] Based on the radiation brightness value, an atmospheric correction algorithm is used in combination with atmospheric parameters during image acquisition, wherein the aerosol optical depth is within a preset aerosol optical depth range in a preset waveband, and the accuracy meets a preset aerosol accuracy standard. The water vapor content is within a preset water vapor content range, and the accuracy meets a preset water vapor accuracy standard. The influence of atmospheric scattering and absorption on the radiation brightness value is eliminated, and the ground reflectance data is obtained, which is within a preset reflectance range, and the accuracy meets a preset reflectance accuracy standard, and the radiation correction is completed.

[0050] Based on the remote sensing image after radiation correction and a high-precision reference map, the reference map adopts a preset coordinate system, and the coordinate accuracy meets a preset map accuracy standard. A ground control point selection algorithm is used to select not less than a preset number of ground control points that are uniformly distributed and have obvious features on the image and the reference map. The pixel coordinates of each ground control point in the image and the geographic coordinates of each ground control point in the reference map are recorded, wherein the pixel coordinate accuracy meets a preset pixel accuracy standard, and the geographic coordinate accuracy meets a preset geographic coordinate accuracy standard.

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

[0052] Based on the coefficient parameters of the polynomial transformation model, the positions of each pixel in the image are recalculated and arranged by a resampling algorithm using a bilinear interpolation method, so that the position deviation of the resampled pixels does not exceed a preset pixel deviation threshold, ensuring that the coordinate system of the image is consistent with the coordinate system of the reference map, and 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 configured to perform initial region identification on the corrected remote sensing image sequence in combination with a pre-trained YOLOv5 model, to obtain an initial region identification range and a corresponding region type of the remote sensing image collected at different time points, and to obtain an edge segmentation feature set corresponding to each initial farmland segmentation sub-region at each time point based on the initial region identification range and an image segmentation algorithm. The region type includes, but is not limited to, farmland, construction land, woodland, pond, and abandoned land.

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

[0055] Based on the corrected remote sensing image sequence, the pre-trained YOLOv5 model is used. The model is trained on a plurality of region sample images including farmland, construction land, woodland, pond, and abandoned land. The corrected single-frame remote sensing image is input into the model. The model's feature extraction network is used to perform multi-scale feature fusion on the image, to obtain feature maps at different levels.

[0056] Based on the feature maps, the prediction head of the YOLOv5 model is used to perform class probability prediction and boundary box regression on each preset anchor box. Anchor boxes with a class probability greater than a preset probability threshold are selected. Anchor boxes with an overlap degree exceeding a preset overlap threshold are merged. Initial identification boundary boxes of each region are obtained. The region type is determined according to the highest class probability corresponding to the anchor box. The initial region identification is completed.

[0057] When the image segmentation unit performs image segmentation and edge feature extraction, based on the initial region identification range and the corresponding region type, a semantic segmentation algorithm is used. The initial region identification range is taken as a region of interest. The image pixels in the region are labeled in detail in terms of class. The farmland is labeled as a pixel set with continuous cultivation texture. The construction land is labeled as a pixel set of artificial building with regular geometric shape. The woodland is labeled as a pixel set with dense vegetation texture. The pond is labeled as a pixel set with uniform water body texture. The abandoned land is labeled as a pixel set without obvious cultivation marks and with sparse vegetation coverage.

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

[0059] It should be further explained that in the embodiment, the gradient change at the boundary of pixels of different categories is calculated, and the pixel points with gradient values exceeding a preset gradient threshold are extracted as edge pixels, and the implementation is specifically as follows:

[0060] Based on the pixel category data after refinement, a preset size square neighborhood window is set for each pixel in the image by an edge detection algorithm. The window includes the center pixel and a preset number of adjacent pixels around the center pixel. Based on the category information of each pixel in the window, the category of the center pixel is compared with the category of each adjacent pixel by a category difference calculation method. The same category is recorded as no difference, and the different category is recorded as difference. The number of adjacent pixels with difference in the window is counted as the initial category difference value of the center pixel.

[0061] Based on the initial category difference value, the gradient components in the horizontal and vertical directions are calculated by a gradient operator. The horizontal direction gradient component is obtained by comparing the category difference between the center pixel and the left and right adjacent pixels. When the category of the left adjacent pixel is different from that of the center pixel, a positive value is assigned, and when the category of the right adjacent pixel is different from that of the center pixel, a negative value is assigned. The sum of the absolute values of the two is taken as the horizontal gradient component. The vertical direction gradient component is obtained by comparing the category difference between the center pixel and the upper and lower adjacent pixels. When the category of the upper adjacent pixel is different from that of the center pixel, a positive value is assigned, and when the category of the lower adjacent pixel is different from that of the center pixel, a negative value is assigned. 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 synthesized by a gradient synthesis algorithm to obtain the gradient value of the center pixel. The gradient value changes positively with the severity of the category 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. The spatial continuity of all candidate edge pixels is detected. By judging whether the spatial distance between adjacent candidate edge pixels is within a preset range, the candidate edge pixels within the preset range are classified into the same edge segment. The edge segments corresponding to the same region type are connected in spatial order to form a closed edge contour. The spatial coordinates and corresponding category difference values of all pixels on the contour are extracted as part of the edge segmentation feature set.

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

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

[0066] It should be further explained that farmland often presents an irregular shape due to topography, farming methods, etc. Traditional center positioning methods are easily affected by shape and lead to deviation. Therefore, the real center is accurately determined by combining multiple algorithms and dynamically adjusting the weight, thereby solving the problem of inaccurate positioning of the center of irregular farmland. At the same time, spatial features and temporal dynamic changes need to be considered in land ecological condition monitoring, and the sub-regions of remote sensing images at different time points lack a unified space-time reference. Therefore, a space-time coordinate system integrating spatial coordinates and time parameters is constructed, and each time point sub-region is embedded therein to realize space-time correlation and unification.

[0067] 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 farmland segmentation sub-region at each time point and a minimum circumscribed rectangle, and obtains a centroid position as a second center point based on the edge pixel position points. The real center position point of the initial farmland segmentation sub-region at each time point is determined based on the mean value of the first center point and the second center point.

[0068] It should be further explained that one implementation of determining the real center position point in the embodiment is as follows:

[0069] Based on the edge pixel space coordinates in the edge segmentation feature set corresponding to the initial farmland segmentation sub-region at each time point, a minimum area circumscribed rectangle algorithm is used to perform multi-directional rotation fitting on all edge pixels, a plurality of rotation directions within a preset angle range are traversed, the area of the rectangle that can contain all edge pixels at each direction is calculated, and the smallest rectangle is selected as the target circumscribed rectangle, which can form a preset angle with the horizontal direction. The four vertex coordinates of the target circumscribed rectangle are extracted, the horizontal coordinate and the vertical coordinate of the intersection point of the two diagonal lines are calculated by the diagonal line midpoint calculation method, and the first center point coordinates are formed by combination, so that the center point is more consistent with the geometric center of the irregular farmland.

[0070] Based on the edge pixel spatial coordinates in the edge segmentation feature set and the category difference information of adjacent pixels, the edge pixels are assigned weight values through a weighted centroid algorithm, wherein the edge pixels that are continuous with the surrounding pixels and have stable category difference are assigned higher weights, and the edge pixels that exist in isolation or have abrupt category difference are assigned lower weights; the horizontal coordinates of all edge pixels are multiplied by the corresponding weights and then summed, and then divided by the sum of all weights to obtain the weighted horizontal average coordinates; the vertical coordinates of all edge pixels are multiplied by the corresponding weights and then summed, and then divided by the sum of all weights to obtain the weighted vertical average coordinates, which are combined to form the coordinates of the second center point, reducing the interference of abnormal pixels in the irregular edge 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 the edge pixels, the average distance of the edge pixels to each side of the minimum area circumscribed rectangle is calculated through a dynamic weight algorithm, and the greater the distance, the more irregular the shape of the farmland, at which time the weight proportion of the coordinates of the second center point is increased, and the smaller the distance, the greater the weight proportion of the first center point; according to the dynamically determined weight, the horizontal coordinates of the first center point and the second center point are weighted and summed according to the weight to obtain the horizontal coordinates of the real center; the vertical coordinates of the two are weighted and summed according to the same weight to obtain the vertical coordinates of the real center, which are combined to form the real center position point of the irregular initial farmland segmentation sub-region at each time point.

[0072] The edge pixels of the irregular farmland are fitted in multiple directions through the minimum area circumscribed rectangle algorithm, and the midpoint of the diagonal of the minimum area circumscribed rectangle is selected as the first center point, which effectively avoids the misjudgment of the geometric center of the irregular region, so that the first center point is more consistent with the actual geometric distribution of the farmland; the edge pixels are assigned differential weights through the weighted centroid algorithm, which reduces the interference of isolated or abrupt edge pixels and improves the representativeness of the second center point coordinates to the edge distribution characteristics; the weight proportions of the two are adaptively adjusted according to the irregularity of the farmland shape by combining the dynamic weight algorithm, so that the finally determined real center position point takes into account both the geometric center and the edge distribution characteristics, significantly improving the accuracy and stability of the center positioning of the irregular farmland, and providing a reliable reference for subsequent feature extraction, spatiotemporal coordinate system construction and other links based on the center position, thereby improving the precision and reliability of the entire land ecological condition monitoring system.

[0073] Based on the real center position point of each time point in the initial farmland segmentation sub-region at the corresponding time, a two-dimensional coordinate is constructed, and the connecting line of all aligned real center position points at consecutive time points is taken as a three-dimensional time axis to construct a spatiotemporal coordinate system, and the initial farmland segmentation sub-region of the remote sensing image collected at each time point is embedded into the spatiotemporal coordinate system at the corresponding time point.

[0074] It should be further pointed out that one construction method of the spatiotemporal coordinate system in the 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, the coordinates of the real center point are subtracted from the spatial coordinates of all edge pixels in the initial cultivated land sub-region using the relative coordinate calculation method, thus obtaining the relative coordinate values ​​of each edge pixel in the local two-dimensional coordinate system and 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 of 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] The process effectively solves the technical problem of inaccurate positioning of the geometric center of irregular farmland by fusing the minimum area circumscribed rectangle algorithm and the weighted centroid algorithm and introducing a dynamic weight adjustment mechanism. The minimum area circumscribed rectangle algorithm accurately captures the macro geometric shape of the farmland by multi-directional rotation fitting, ensuring that the first center point conforms to the actual distribution. The weighted centroid algorithm suppresses the interference of abnormal edge pixels by differentiating the weights, improving the representation ability of the coordinates of the second center point to edge details. Further, the weights of the two are adaptively distributed according to the irregularity of the shape, so that the finally determined real center position point has both geometric inclusiveness and distribution representativeness, significantly improving the accuracy and robustness of the center positioning. On this basis, by constructing a three-dimensional space-time coordinate system integrating spatial coordinates and time dimension, the multi-temporal farmland sub-regions are uniformly embedded into the space-time framework, realizing the time sequencing and spatial integration management of the monitoring data. The space-time coordinate system ensures the spatial comparability and temporal continuity of multi-temporal data through local two-dimensional coordinates combined with three-dimensional time axis and plane vertical adaptation, providing a stable and reliable space-time reference for subsequent analysis such as feature extraction and change detection. Finally, the technical scheme enhances the perception ability and monitoring accuracy of irregular farmland changes by improving the center positioning accuracy and constructing a unified space-time reference.

[0084] The feature extraction unit extracts features along any N edge extension directions in the initial farmland segmentation sub-region corresponding to each time point, taking the real center position point of the initial farmland segmentation sub-region corresponding to the time point as the starting point, to obtain a spatial feature sequence and a spatial feature fluctuation value of any N edge extension directions at the current time point.

[0085] It should be further noted that in the embodiment, one implementation of the feature extraction along any N edge extension directions in the initial farmland segmentation sub-region corresponding to the time point is specifically as follows:

[0086] Based on the edge segmentation feature set of the initial farmland segmentation sub-region at each time point, the edge segmentation feature set includes the spatial coordinates of the edge pixels and the category difference information of the adjacent pixels. The gradient direction of each edge pixel is quantized by a direction clustering algorithm, wherein the gradient direction is the direction of the current pixel pointing to the adjacent but different category pixel.

[0087] The gradient direction distribution of all edge pixels is counted, and the edge pixels with a direction difference within a preset range are classified into the same direction cluster. The number of edge pixels included in each direction cluster is calculated.

[0088] Selecting a preset number of direction clusters with the largest number of pixels, taking the average of the gradient directions of each direction cluster as the representative direction of the cluster, and performing angle calibration on the representative directions so that the included angle difference between any two adjacent representative directions does not exceed a preset angle deviation, to obtain the angle parameters of any N edge extension directions;

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

[0090] By a pixel traversal algorithm, each pixel's spatial coordinates and the attribute of the sub-region to which the pixel belongs are recorded in sequence along the ray direction according to the pixel arrangement order, and the pixels contacted by the ray after the ray exits the initial farmland segmentation sub-region are removed, forming a pixel sequence arranged continuously from the center to the edge in each direction, and the sequence length is the total number of pixels penetrated 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, wherein the normalized vegetation index is obtained by the operation of the near-infrared band reflectance and the red band reflectance, and the enhanced vegetation index is obtained by the operation of the near-infrared band reflectance, the red band reflectance and the blue band reflectance; the normalized vegetation index value and the enhanced vegetation index value of each pixel are recorded in sequence according to the arrangement order of the pixel in the sequence, forming the normalized vegetation index sequence and the enhanced vegetation index sequence of each edge extension direction, and the numerical position in the sequence corresponds to the position of the corresponding pixel in the ray direction;

[0092] Based on the normalized vegetation index sequence and the enhanced vegetation index sequence in each direction, the normalized vegetation index values of two adjacent pixels in the sequence are differentiated by a difference calculation algorithm, and the absolute value of the operation result is taken as the normalized vegetation index spatial change value of the adjacent pixel pair;

[0093] The enhanced vegetation index values of two adjacent pixels in the sequence are differentiated in the same way, and the absolute value of the operation result is taken as the enhanced vegetation index spatial change value of the adjacent pixel pair;

[0094] According to the appearance order of the adjacent pixel pairs in the pixel sequence, all the spatial change values are arranged in sequence to form a spatial feature fluctuation value sequence adapted to the length of the spatial feature sequence, each fluctuation value corresponds to the difference between two adjacent pixels in the spatial feature sequence, and finally the spatial feature sequence contained in the feature space of the remote sensing image at the current time point and the spatial feature fluctuation value are obtained.

[0095] The process significantly improves the comprehensiveness and accuracy of capturing the spatial characteristics of cultivated land vegetation by combining multi-directional ray scanning with gradient direction clustering feature extraction method. The edge pixel gradient direction clustering analysis autonomously determines the most representative edge extension direction, overcoming the subjectivity and limitations of manual selection of direction, ensuring that the feature extraction direction is highly consistent with the actual contour characteristics of cultivated land. By emitting rays from the true center position point to each edge direction and traversing the pixel sequence, the gradual feature collection from the core to the edge is realized, effectively capturing the spatial heterogeneity and transition characteristics 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 of vegetation density and activity, but also accurately quantify the boundary mutation and internal gradual change information of vegetation characteristics through the spatial feature fluctuation value generated by the difference operation of adjacent pixels. This multi-directional, multi-index, and multi-scale feature extraction system can comprehensively represent the spatial distribution pattern and change trend of cultivated land vegetation, providing high-precision and high-reliability feature data support for subsequent cultivated land area change monitoring and ecological degradation inversion, greatly enhancing the perception ability and early warning accuracy of the system for complex-shaped cultivated land and gradual ecological change.

[0096] The region determination module determines the true cultivated land segmentation sub-region area based on the remote sensing image feature space combined with the initial cultivated land segmentation sub-region, and determines the cultivated land area fluctuation value based on the determined true cultivated land segmentation sub-region area change value; it needs to be further explained that the region determination module in the embodiment includes a region determination unit, a region area unit, and a region fluctuation determination unit;

[0097] The region determination unit determines the second edge position point set based on the spatial feature sequence and the spatial feature fluctuation value of any N edge extension directions at the current time point, and the region type fluctuation threshold value. When the spatial feature fluctuation value of any edge extension direction at M consecutive pixel position points does not satisfy the region type fluctuation threshold value, the pixel position point corresponding to the maximum spatial feature fluctuation value is taken as the edge point of the current edge extension direction. Based on this, the N times are repeated to determine the second edge position point set.

[0098] It needs to be further explained that the specific process of determining the second edge position point set in the embodiment includes:

[0099] Based on the historical remote sensing images of the same growth season in the target region for consecutive years, the images are segmented into sub-regions such as cultivated land, construction land, forest land, pond, and abandoned land through a region type classification algorithm, and the spatial feature fluctuation values of all edge extension directions in each sub-region are extracted to form a type-specific historical sample set;

[0100] For the sample set of the same region type, the statistical analysis algorithm is used to calculate the arithmetic mean of the spatial feature fluctuation value, and the deviation degree of all values in the sample set from the mean value is calculated to obtain a statistical quantity reflecting the dispersion characteristic.

[0101] Based on the mean and the dispersion statistics, and combined with the inherent properties of the region type, the threshold calibration algorithm is used to form the spatial feature fluctuation threshold of the region type by subtracting the preset multiple of the dispersion statistics from the mean as the lower limit and adding the preset multiple of the dispersion statistics to the mean as the upper limit, wherein the lower limit of the fluctuation threshold of the cultivated land type needs to be higher than the typical fluctuation value of the abandoned land, and the upper limit needs to be lower than the sudden fluctuation value of the construction land; for example, the inherent properties of the region type, such as the seasonal fluctuation of the cultivated land due to crop growth, the fluctuation threshold needs to cover the normal change range from the green-up period to the mature period; the fluctuation threshold of the construction land is set to a smaller range to limit abnormal changes due to the stable structure;

[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 true center position point to the edge) and the fluctuation threshold of the corresponding cultivated land type, the point-by-point traversal algorithm is used to start from the first pixel position point in the sequence and check in turn, 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, mark it as a pixel that does not meet the fluctuation threshold, and at the same time start the continuous counting.

[0103] If the subsequent adjacent pixels are continuously pixels that do not meet the fluctuation threshold, the count is incremented until the continuous 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 is stopped, the M pixel position points are arranged in a continuous sub-sequence according to the arrangement order, and the start position and the end position of the sub-sequence are recorded.

[0104] Based on the spatial feature fluctuation values of the locked continuous M pixel position points, the point-by-point comparison algorithm is used to compare the fluctuation values of each pixel in the sub-sequence in turn, and filter out the pixel position point with the maximum fluctuation value.

[0105] If there are multiple pixels with the same fluctuation value and all of them are the maximum value, the pixel position point closest to the edge direction in the sub-sequence is selected as the target point; the target point is marked as the edge point of the current edge extension direction, because its maximum fluctuation value corresponds to the critical mutation position from the internal feature to the external feature, the boundary of the cultivated land sub-region can be accurately defined.

[0106] The above detection, locking and screening operations are performed on any N edge extension directions respectively, and an edge point is obtained for each extension direction; through the spatial coordinate verification algorithm, it is checked whether the spatial distribution of all edge points conforms to the contour trend of the initial cultivated land segmentation sub-region, and the abnormal edge points deviating from the overall trend are removed, if there is no edge point meeting the conditions in a certain direction, the coordinates of the edge points in the adjacent directions are interpolated to supplement; all the edge points after verification are connected in turn according to the angle order of the edge extension directions to form a closed edge contour, and the spatial coordinates of all the edge points included are summarized as the second edge position point set.

[0107] The process effectively improves the accuracy and adaptability of regional edge detection by comprehensively utilizing historical remote sensing image data and spatial feature fluctuation analysis. The core lies in constructing a fluctuation threshold range of regional types based on multi-temporal images, fully considering the change characteristics of different ground object properties, so that the edge judgment has both statistical significance and conforms to the ground object evolution law. By detecting the pixel sequence that continuously does not satisfy the fluctuation threshold in each direction and locating the fluctuation extreme points, the feature mutation position can be accurately captured, the misjudgment caused by noise or local change can be significantly suppressed, and the position of the edge point is ensured to be in the real transition zone between the internal feature and the external environment. Further, with the help of spatial verification and interpolation mechanism, the spatial consistency and contour closure of the edge point are ensured, and the final second edge position point set not only has clear boundaries and complete structure, but also has high land class distinguishing ability, especially suitable for fine definition of farmland and abandoned land, construction land and other easily confused areas, providing reliable technical support for remote sensing image segmentation and land use dynamic monitoring.

[0108] Based on the second edge position point set of the initial farmland segmentation sub-region at each time point and the corresponding initial farmland segmentation sub-region, the regional integral algorithm is used to obtain the area of the real farmland segmentation sub-region corresponding to each time point, and the pixel point corresponding to the largest real farmland segmentation sub-region is taken as the starting point, and the absolute value of the time feature change value of the adjacent real farmland segmentation sub-region area at the same pixel position point is extracted along the positive and negative axes of the three-dimensional time axis.

[0109] It should be further pointed out that the specific process of obtaining the area of the real farmland segmentation sub-region in the embodiment includes:

[0110] Based on the second edge position point set of the initial farmland segmentation sub-region at each time point and the corresponding initial farmland segmentation sub-region, the boundary verification algorithm is used to compare the spatial position of each edge point in the second edge position point set with the edge pixels of the initial farmland segmentation sub-region, eliminate the abnormal edge points beyond the range of the initial sub-region, and retain the effective edge points within the range of the initial sub-region; the effective edge points are smoothed according to the spatial distribution to eliminate the abrupt jumps of adjacent points, and a continuous boundary transition point set is formed;

[0111] Based on the smoothed effective edge points, the polygon fitting algorithm is used to connect all the effective edge points in the clockwise or counterclockwise order of the edge extension direction, form a closed boundary contour surrounding the farmland segmentation sub-region, and ensure 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 area within the closed boundary contour is divided into grid units consistent with the pixel size of the image by 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 the length.

[0113] When the area unit performs area integration, the area integration algorithm is used to count all grid units within the closed boundary contour. The number of grid units completely contained in the contour is counted. For grid units partially contained in the contour, the proportion of the area falling within the contour is converted. For example, if the falling proportion exceeds a preset threshold, it is counted as a complete unit, otherwise it is counted as a partial unit according to the proportion.

[0114] The areas of all completely contained grid units and the converted areas of partial units are summed up to obtain the real farmland segmentation sub-area area corresponding to each time point.

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

[0116] Based on the real farmland segmentation sub-area area at consecutive time points, the maximum area real farmland segmentation sub-area is selected as the reference sub-area by the extreme value screening algorithm by comparing the areas of all time points. The spatial coordinates of all pixel position points in the reference sub-area are recorded to form a reference pixel point set.

[0117] Based on the time parameter sorting of the three-dimensional time axis, the time point corresponding to the reference sub-area is taken as the origin. The direction definition rule is used to set the direction extending from the origin to the past time point as the negative axis direction of the three-dimensional time axis, and the direction extending to the future time point as the positive axis direction, to clearly arrange the continuous time points on the positive and negative axes.

[0118] Based on each pixel position point in the reference pixel point set, the coordinate mapping algorithm is used to map its spatial coordinates to the real farmland segmentation sub-area corresponding to the adjacent time points in the positive and negative axes of the three-dimensional time axis, respectively. The pixel position points in each adjacent sub-area with the same spatial coordinates as those in the reference pixel point set are found to form a corresponding relationship set of the same pixel position at different time points. The pixel points that cannot be matched to the same coordinates are marked as missing pixel points.

[0119] Based on the positive and negative axis directions of the three-dimensional time axis, a sequence traversal algorithm is used to select, from a reference sub-region corresponding time point, in turn, a previous sequence adjacent time point along the negative axis direction and a subsequent sequence adjacent time point along the positive axis direction, to form a plurality of groups of adjacent time point pairs, including the reference time point and the previous sequence time point, the reference time point and the subsequent sequence time point, adjacent previous sequence time points, and adjacent subsequent sequence time points;

[0120] Based on the corresponding relationship of the same pixel position in the adjacent time point pair, a difference value operation algorithm is used to extract the area proportion of the pixel position in the real cultivated land segmentation sub-region corresponding to the adjacent two time points, that is, the area proportion contributed by the pixel point in the sub-region, and the area proportion difference value is obtained by subtracting the area proportion of the previous time point from the area proportion of the subsequent time point.

[0121] Based on the absolute value conversion algorithm, the area proportion difference value is processed, and the absolute value thereof is taken as the absolute value of the time feature change value of the pixel position in the corresponding adjacent time point pair. For the position marked as a missing pixel point, a preset missing value is assigned as the absolute value of the time feature change value, and finally an absolute value set of the time feature change value of all the same pixel positions at continuous time points is formed.

[0122] The regional fluctuation determination unit is configured to obtain a first cultivated land area fluctuation value sequence according to the real cultivated land segmentation sub-region areas corresponding to adjacent time points, and obtain the proportion of the absolute value of the time feature change value of the same pixel position that does not satisfy the regional type fluctuation threshold as a second cultivated land area fluctuation value sequence based on the absolute value of the time feature change value and the regional type fluctuation threshold. The cultivated land area fluctuation value of the adjacent time points is obtained based on the mean value of the first cultivated land area fluctuation value and the second cultivated land area fluctuation value of the corresponding adjacent time points in the first cultivated land area fluctuation value sequence and the second cultivated land area fluctuation value sequence.

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

[0124] Based on the real cultivated land segmentation sub-region areas corresponding to continuous time points and adjacent time point pairs (each time point and the next time point are grouped into an adjacent time point pair according to the time sequence on the three-dimensional time axis), a difference value operation method is used to calculate the difference value between the real cultivated land segmentation sub-region area of the subsequent time point and the real cultivated land segmentation sub-region area of the previous time point in each adjacent time point pair, and the absolute value of the difference value 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 turn to form a first cultivated land area fluctuation value sequence, and each fluctuation value in the sequence corresponds to a corresponding adjacent time point pair.

[0125] Based on the absolute value set of the time feature change value of the same pixel position point at continuous time points (obtained by the area unit along the positive and negative axes of the three-dimensional time axis) and the regional type fluctuation threshold corresponding to the cultivated land type (obtained by statistical calibration of the spatial feature fluctuation value of the historical cultivated land sample), the absolute value of the time feature change value of all pixel position points in each adjacent time point pair is checked one by one through the threshold comparison algorithm. If the absolute value of the time feature change value of a pixel position point is greater than the upper limit of the regional type fluctuation threshold or less than the lower limit, the pixel position point is determined to be a pixel position point that does not satisfy the regional type fluctuation threshold. The number of pixel position points that do not satisfy the regional type fluctuation threshold in the adjacent time point pair is calculated by the proportion statistical algorithm, and the ratio of the number to the total number of pixel position points in the adjacent time point pair is taken as the second cultivated land area fluctuation value of the adjacent time point pair. The total pixel position points are all pixel points in the reference sub-region in the adjacent sub-region, and do not include missing pixel points. According to the arrangement order of the adjacent time point pair on the three-dimensional time axis, all second cultivated land area fluctuation values are arranged in turn to form a second cultivated land area fluctuation value sequence, and each fluctuation value in the sequence corresponds to an adjacent time point pair.

[0126] Based on the first cultivated land area fluctuation value sequence and the second cultivated land area fluctuation value sequence (both sequences have the same length and are consistent with the number of adjacent time point pairs), the first cultivated land area fluctuation value and the second cultivated land area fluctuation value corresponding to the same adjacent time point pair in the two sequences are summed by the mean value calculation method, and the sum is divided by two to obtain the average fluctuation value of the adjacent time point pair. The average fluctuation value is taken as the cultivated land area fluctuation value of the adjacent time point. The mean value calculation process is repeated for all adjacent time point pairs to obtain a set of cultivated land area fluctuation values of continuous adjacent time points, and each fluctuation value in the set corresponds to a group of adjacent time point pairs on the three-dimensional time axis.

[0127] The process effectively realizes high-precision dynamic monitoring of cultivated land area and its volatility by fusing multi-temporal spatial boundary information and time series change characteristics. Its advantage lies in the use of a double-checking mechanism of real edge point set and initial region, which significantly improves the accuracy of area calculation and ensures that the integral result is more in line with the actual distribution of ground objects. By introducing pixel-level change tracking in the positive and negative directions of the three-dimensional time axis, it can sensitively capture the area ratio difference at the same location at different time points, thereby extracting time characteristic change values with physical meaning. Further combining macro-area fluctuations with micro-pixel changes, an integrated fluctuation index is constructed through mean fusion, reflecting both the overall area change of the region and the abnormal change information of the internal pixels, making the final obtained cultivated land area fluctuation value both robust and sensitive. This method effectively suppresses false positives caused by temporary interference or local noise and can accurately identify sustained changes caused by real farming activities, providing reliable data for dynamic monitoring of cultivated land resources, abandoned land warning, and land use assessment.

[0128] The inversion module obtains the ecological degradation fluctuation value of the real cultivated land segmentation sub-region based on the feature space of remote sensing images combined with the land surface temperature and the inversion of the cultivated land vegetation coverage.

[0129] It should be further explained that the inversion module in the embodiment includes a coverage inversion unit, a first change trend unit, a second change trend unit, and a comprehensive trend unit.

[0130] The coverage inversion unit inversely obtains the corresponding cultivated land vegetation coverage and the corresponding coverage distribution probability function of each real cultivated land segmentation sub-region at the corresponding time point based on the normalized vegetation index sequence and the enhanced vegetation index sequence corresponding to each time point and the normalized vegetation index spatial change value and the enhanced vegetation index spatial change value at the corresponding time point. Meanwhile, the coverage inversion unit obtains the time distribution function of the cultivated land vegetation coverage based on the corresponding cultivated land vegetation coverage of each real cultivated land segmentation sub-region and the absolute value of the time characteristic change value of the area of the adjacent real cultivated land segmentation sub-region at the same pixel position at consecutive time points.

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

[0132] Based on the known sample area of cultivated land vegetation coverage, the sample area needs to include cultivated land plots with high, medium, and low coverage levels and have field-measured cultivated land vegetation coverage data. Through a sample calibration algorithm, the normalized vegetation index value, enhanced vegetation index value, and corresponding measured coverage value of each pixel in the sample area are extracted 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 sequence, the enhanced vegetation index sequence and the global nonlinear mapping model at each time point, the initial coverage value of each pixel in the sub-region of the real cultivated land is obtained by the pixel-by-pixel analysis algorithm, the corresponding values of the pixel in the two index sequences are extracted, and the global nonlinear mapping model is input synchronously; the initial coverage values of adjacent pixels are compared by the spatial consistency verification algorithm, if the difference exceeds the preset consistency threshold, the coverage values of surrounding pixels are interpolated and corrected to ensure the continuity of spatial distribution; the corrected coverage values of all pixels are arranged according to their spatial coordinates to form the cultivated land vegetation coverage which is a basic unit of pixels and presents spatial heterogeneity.

[0138] Based on the cultivated land vegetation coverage, the coverage values are divided into several continuous intervals according to the preset interval by the interval division algorithm, wherein the interval size is determined according to the standard deviation of the sample coverage; the proportion of the weighted pixel number in the total weighted pixel number is taken as the interval probability value by the weighted statistical algorithm, according to the distance of the pixel in each interval from the center of the sub-region, the closer the distance, the higher the weight, which reflects the representativeness of the core area; the coverage distribution probability function is obtained by fitting the maximum likelihood estimation method with the coverage interval as the horizontal axis and the interval probability value as the vertical axis by the function optimization algorithm, and the function needs to be verified by the chi-square test to ensure that the fitting degree can accurately reflect the spatial probability distribution characteristics of the coverage.

[0139] Based on the cultivated land vegetation coverage at continuous time points and the absolute value of the time characteristic change value of the same pixel position point, all time points are sorted according to the collection time by the time axis alignment algorithm to ensure that the interval of adjacent time points is uniform.

[0140] The coverage time sequence of each pixel position point is formed by the pixel-level time correlation algorithm, the coverage values of the pixel at continuous time points are extracted, and the absolute value of the corresponding time characteristic change value is taken as the weight of each time point in the sequence (when the value is zero, a preset basic weight is taken); the coverage time sub-function of each pixel is obtained by the dynamic trend fitting algorithm, the weighted least squares method is used to fit the trend of the coverage time sequence of each pixel; the coverage time sub-function of all pixels is spliced according to its spatial coordinates by the spatial integration algorithm to form a binary function containing spatial position parameters and time parameters, that is, the cultivated land vegetation coverage time distribution function, which can output the coverage prediction value of any time point and any spatial position.

[0141] The first change trend unit is configured to, based on the corresponding farmland vegetation coverage and the corresponding coverage distribution probability function of each real farmland segmentation sub-region, combine the ground surface temperature and a thermal inertia model to inversely derive a moisture content distribution state function of the corresponding real farmland segmentation sub-region, and based on the moisture content distribution state function of the real farmland segmentation sub-region, combine a farmland vegetation coverage time distribution function to obtain a moisture content time fluctuation function corresponding to the real farmland segmentation sub-region at a continuous time point;

[0142] It should be further explained that one implementation of the moisture content distribution state function of the real farmland segmentation sub-region in the embodiment includes:

[0143] Based on the remote sensing image at the corresponding time point, the thermal radiation value of each pixel in the real farmland segmentation sub-region is extracted through a thermal infrared band analysis algorithm, the thermal radiation value is converted into a ground surface temperature value, and ground surface temperature distribution data containing the spatial coordinates of all pixels and the corresponding temperature values are formed. At the same time, based on the farmland vegetation coverage, the ground surface temperature distribution data is corrected for vegetation coverage, and a vegetation shading coefficient is assigned to the temperature value of a high-coverage pixel, wherein the higher the coverage, the larger the vegetation shading coefficient, which reduces the interference of vegetation on soil temperature, and corrected ground surface temperature distribution data is obtained.

[0144] Based on the corrected ground surface temperature distribution data and a thermal inertia model, wherein the model reflects the correlation between soil thermal inertia and temperature daily range, thermal conduction characteristics, the larger the thermal inertia, the slower the soil temperature changes with the environment, and through a day and night temperature difference algorithm, the ground surface temperature daily range of each pixel in the sub-region is obtained, specifically the difference between the daytime maximum temperature and the nighttime minimum temperature. Combined with the shortwave albedo data of the remote sensing image, the thermal inertia model is used to substitute the temperature daily range and the albedo into the calculation to obtain the thermal inertia value of each pixel, and thermal inertia spatial distribution data is formed.

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

[0146] Based on the coverage distribution probability function and the vegetation coverage of cultivated land, the sub-regions are divided into several layers, including high coverage layer, medium coverage layer and low coverage layer, by hierarchical correction algorithm according to the coverage interval, and 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, the vegetation transpiration correction coefficient is given to the high coverage layer to reduce the sensitivity of thermal inertia to water content, and the soil bare correction coefficient is given to the low coverage layer to improve the sensitivity of thermal inertia to water content, and the local mapping relationship of each layer is obtained.

[0147] Based on the thermal inertia spatial distribution data, the local mapping relationship of each layer and the coverage distribution probability function, the initial water content value of each pixel is converted from the thermal inertia value by spatial interpolation algorithm according to the coverage value of the pixel to determine the corresponding layer level, and the local mapping relationship of the corresponding layer level is called; the initial water content value is adjusted by probability weighting algorithm combined with the probability value corresponding to the coverage of the pixel in the coverage distribution probability function, wherein the higher the probability, the smaller the adjustment range, and the result is ensured to meet the overall distribution characteristics.

[0148] The adjusted water content values of all pixels are integrated according to the spatial coordinates to form a water content distribution state function describing the spatial distribution law of water content in the sub-region, and the function input is the spatial coordinates of the pixel and the output is the corresponding water content prediction value.

[0149] It should be further pointed out that one implementation of the water content time fluctuation function corresponding to the real cultivated land segmented sub-region at continuous time points in the embodiment includes:

[0150] Based on the acquisition time of continuous time points and the spatial coordinates of the real cultivated land segmented sub-region, all time points are arranged in chronological order by time axis calibration algorithm to ensure the uniformity of the interval of adjacent time points, and if there is interval difference, interpolation is completed to form an equidistant time sequence axis, so that the water content data at different time points have time comparability.

[0151] Based on the spatial coordinates of each pixel in the real cultivated land segmented sub-region and the water content distribution state function at each time point, the water content value of the pixel at each time point is extracted by space-time indexing algorithm, arranged in order of the time sequence axis to form the water content time sequence of the pixel; for the missing water content values in the time sequence, the adjacent time point water content values and the prediction values of the cultivated land vegetation coverage time distribution function of the pixel are interpolated to fill in to ensure the continuity of the sequence.

[0152] Based on the time distribution function of the cultivated land vegetation coverage, the time change rate of the coverage of each pixel at continuous time points is extracted, reflecting the increase and decrease amplitude of the coverage over time. Through a conversion algorithm, the time change rate of the coverage is converted into the weight value of the time sequence of the water content, wherein the greater the time change rate of the coverage, the higher the weight, reflecting the significant influence of vegetation change on the fluctuation of the water content.

[0153] Based on the time sequence of the water content of each pixel and the corresponding weight value, through a weighted time sequence analysis algorithm, the trend decomposition is performed on the water content values in the sequence, and the long-term trend component and the short-term fluctuation component are separated.

[0154] Through a function fitting algorithm, the time is taken as the independent variable, the water content value is taken as the dependent variable, the periodic characteristics of the short-term fluctuation component are combined, such as the daily change and the weekly change rule, and the weighted least square method is used to fit to obtain the water content time fluctuation sub-function of the pixel.

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

[0156] Based on the cultivated land vegetation coverage corresponding to each real cultivated land segmented sub-region and the corresponding coverage distribution probability function, and in combination with the water content distribution state function of each real cultivated land segmented sub-region, the soil fertility distribution function of each real cultivated land segmented sub-region is obtained through inversion. According to the soil fertility distribution function of each real cultivated land segmented sub-region at continuous time points and the cultivated land vegetation coverage time distribution function, the first soil fertility time change trend at continuous time points is obtained through inversion.

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

[0158] Firstly, based on the spatial change characteristics of the normalized vegetation index sequence and the enhanced vegetation index sequence, through a vegetation growth trend extraction algorithm, the difference value between the normalized vegetation index value and the enhanced vegetation index value of each pixel at adjacent time points is calculated, and the change rate value is obtained by dividing the difference value by the time interval; according to the threshold range determined in advance through historical data statistics, the change rate value is divided into three discrete levels, i.e., a high change rate level, a medium change rate level and a low change rate level; based on the priori knowledge of agricultural ecology, a corresponding rule between the level and the fertility is established, that is, the high change rate level represents vigorous vegetation growth, corresponding to a high soil basic fertility level, the medium change rate level represents stable vegetation growth, corresponding to a medium soil basic fertility level, and the low change rate level represents slow vegetation growth, corresponding to a low soil basic 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] The slope value of the coverage sequence is calculated as a long-term change characteristic by linear regression analysis combined with the trend characteristics of the time distribution function of the cultivated land vegetation coverage. A change amount correction algorithm is adopted. When the slope is negative, the initial change amount is multiplied by a decay coefficient less than 1. When the slope is positive, the initial change amount is multiplied by an enhancement coefficient greater than 1. When the slope is close to zero, the initial change amount is kept unchanged. Finally, the change amounts of all pixels are arranged in time sequence by a spatial aggregation algorithm, and the average value of the change amounts of all pixels at each time point is calculated to form the first soil fertility time change trend representing the overall change trend of the regional fertility.

[0164] The second change trend unit obtains the crop phenological node corresponding to the real cultivated land segmentation sub-region based on the normalized vegetation index spatial change value and the enhanced vegetation index spatial change value combined with the crop phenological period extraction algorithm, performs phenological period anomaly identification based on the time characteristic change value corresponding to the crop phenological node corresponding to the real cultivated land segmentation sub-region, obtains the corresponding crop phenological period anomaly identification feature space, and inversely obtains the second soil fertility time change trend of the corresponding real cultivated land segmentation sub-region based on the corresponding crop phenological period anomaly identification feature space.

[0165] It should be further explained that the process of obtaining the crop phenological node corresponding to the real cultivated land segmentation sub-region in the embodiment includes:

[0166] Based on the normalized vegetation index spatial change value and the enhanced vegetation index spatial change value, the normalized vegetation index spatial change value of each pixel at consecutive time points is extracted by a time sequence reorganization algorithm, and the normalized vegetation index spatial change time sequence of the pixel is formed by arranging the values in time sequence. Meanwhile, the enhanced vegetation index spatial change value of each pixel at consecutive time points is extracted, and the enhanced vegetation index spatial change time sequence of the pixel is formed by arranging the values in time sequence. The normalized vegetation index spatial change time sequences of all pixels in the real cultivated land segmentation sub-region are aligned by a regional aggregation algorithm, and the arithmetic average value of the corresponding values of all pixels at each time point is calculated to obtain the normalized vegetation index spatial change time sequence at the sub-region level. The same method is used to obtain the enhanced vegetation index spatial change time sequence at the sub-region level.

[0167] Based on the two types of sub-region level index spatial variation time series, through the sliding average filtering algorithm, a fixed length of three time points window is adopted, for each non-end point position of the time point in the sequence, the arithmetic mean value of the numerical value of the previous time point, the current time point and the next time point is calculated, and the average value is assigned to the current time point as the smoothing value; the original value of the time point at both ends of the sequence is retained; through the trend consistency verification algorithm, the linear regression slope difference value of the sequence before and after smoothing in the overall time range is calculated, if the absolute value of the slope difference value exceeds the preset threshold, the window length is reduced by one time point and the sliding average calculation is performed again, until the slope difference value meets the requirements, and finally the denoised normalized vegetation index spatial variation smoothing sequence and the enhanced vegetation index spatial variation smoothing sequence are obtained.

[0168] Based on the two types of sub-region level index spatial variation time series, through the sliding average filtering algorithm, a fixed length of three time points window is adopted, for each non-end point position of the time point in the sequence, the arithmetic mean value of the numerical value of the previous time point, the current time point and the next time point is calculated, and the average value is assigned to the current time point as the smoothing value; the original value of the time point at both ends of the sequence is retained; through the trend consistency verification algorithm, the linear regression slope difference value of the sequence before and after smoothing in the overall time range is calculated, if the absolute value of the slope difference value exceeds the preset threshold, the window length is reduced by one time point and the sliding average calculation is performed again, until the slope difference value meets the requirements, and finally the denoised normalized vegetation index spatial variation smoothing sequence and the enhanced vegetation index spatial variation smoothing sequence are obtained.

[0169] Based on the two types of sub-region level index spatial variation time series, through the sliding average filtering algorithm, a fixed length of three time points window is adopted, for each non-end point position of the time point in the sequence, the arithmetic mean value of the numerical value of the previous time point, the current time point and the next time point is calculated, and the average value is assigned to the current time point as the smoothing value; the original value of the time point at both ends of the sequence is retained; through the trend consistency verification algorithm, the linear regression slope difference value of the sequence before and after smoothing in the overall time range is calculated, if the absolute value of the slope difference value exceeds the preset threshold, the window length is reduced by one time point and the sliding average calculation is performed again, until the slope difference value meets the requirements, and finally the denoised normalized vegetation index spatial variation smoothing sequence and the enhanced vegetation index spatial variation smoothing sequence are obtained.

[0170] Based on the two types of sub-region level index spatial variation time series, through the sliding average filtering algorithm, a fixed length of three time points window is adopted, for each non-end point position of the time point in the sequence, the arithmetic mean value of the numerical value of the previous time point, the current time point and the next time point is calculated, and the average value is assigned to the current time point as the smoothing value; the original value of the time point at both ends of the sequence is retained; through the trend consistency verification algorithm, the linear regression slope difference value of the sequence before and after smoothing in the overall time range is calculated, if the absolute value of the slope difference value exceeds the preset threshold, the window length is reduced by one time point and the sliding average calculation is performed again, until the slope difference value meets the requirements, and finally the denoised normalized vegetation index spatial variation smoothing sequence and the enhanced vegetation index spatial variation smoothing sequence are obtained.

[0171] It needs to be further explained that the process of obtaining the second soil fertility time change trend corresponding to the segmented sub-region of the real cultivated land in the embodiment includes:

[0172] Based on the crop phenological nodes corresponding to the segmented sub-region of the real cultivated land at the historical time points of more than five consecutive years, the arithmetic mean of the dates of the green-up stage, the jointing stage, the grain-filling stage, and the maturation stage nodes in the same agricultural year in previous years is calculated as the average time value by time series statistical algorithm, and the standard deviation of these date values is calculated as the time fluctuation range. Three key feature parameters are extracted from the historical normalized difference vegetation index (NDVI) spatial change value sequence at the same period by multi-feature extraction algorithm: the peak value, i.e., the maximum value of the sequence, the peak duration, i.e., the number of consecutive time points at which the value remains above 80% of the peak value, and the curve rising rate, i.e., the average daily change amount from the growth starting point to the peak point. Meanwhile, the same three feature parameters are extracted from the enhanced vegetation index (EVI) spatial change value sequence, forming a dual-phenological-period index feature benchmark set containing time features and index features.

[0173] The specific dates of the current phenological nodes are converted into annual day serial numbers, i.e., the day of the year, by date conversion algorithm. The average time day serial number of the corresponding phenological nodes in the benchmark set is subtracted from the current day serial number by difference calculation algorithm, obtaining the original time deviation value. The time feature change value of each pixel position is converted into a weight coefficient in the range of 0 to 1 using the min-max normalization method by weight calculation algorithm, and the weight coefficient is multiplied by the original time deviation value to obtain the weighted phenological 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. A deviation value exceeding the upper limit is marked as a delay anomaly, and a deviation value below the lower limit is marked as an advance anomaly. The relative deviation percentage of 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 by feature comparison algorithm. When the relative deviation percentage is below -30%, it is marked as a lack of vitality anomaly. All anomaly markers are integrated by spatial aggregation algorithm according to the type of phenological period and pixel coordinates to form a crop phenological period anomaly recognition feature space containing anomaly type, geographic coordinates, and quantitative deviation degree.

[0175] Based on the crop phenology anomaly identification feature space and historical soil fertility change sample data, a quantitative relationship model between the phenology anomaly features and the soil fertility change amplitude is established by a multiple linear regression algorithm: the delay anomaly days of the green-up period are taken as continuous independent variable one, the product of the delay anomaly days of the jointing period and the insufficient vigor anomaly degree is taken as continuous independent variable two, the weighted sum of the peak anomaly degree of the grain-filling period and the delay anomaly days of the maturity period is taken as continuous independent variable three, and the soil fertility change amplitude is taken as continuous dependent variable, and the regression coefficients are solved by the least square method; the regression coefficients are assigned with weight coefficients according to the absolute value size for different anomaly types to form a weighted phenology-fertility correlation model.

[0176] The anomaly types and deviation degree values in the crop phenology anomaly identification feature space are input into the weighted phenology-fertility correlation model by a change amount calculation algorithm to calculate the soil fertility change amount corresponding to each phenology node; the fertility change amounts of the nodes are connected in the order of the phenology time to form an initial change sequence by a time sequence construction algorithm; the initial change sequence is processed by a moving average with a window size of three time points to eliminate short-term fluctuations and extract long-term trend components by a time sequence smoothing algorithm; the slope values of the spatial change trends of the normalized vegetation index and the enhanced vegetation index are calculated by a trend verification algorithm, and when the slope of the normalized vegetation index and the slope of the enhanced vegetation index are opposite in sign, the fertility change trend is adjusted by weighting, and the adjustment weight is the absolute value of the ratio of the two slopes, and finally the second soil fertility time change trend of the real farmland segmentation sub-region is formed.

[0177] The trend unit is used for time alignment processing of the first soil fertility time change trend and the second soil fertility time change trend, and the ecological degradation fluctuation value of the real farmland segmentation sub-region is obtained by a combined weighted average algorithm of the first soil fertility time change trend and the second soil fertility time change trend in the aligned time period.

[0178] The process builds a comprehensive monitoring system from vegetation dynamics, soil moisture to soil fertility changes by deeply fusing multi-source remote sensing data and crop ecological parameter data, significantly improving the comprehensiveness and accuracy of cultivated land ecological state evaluation. Its core advantage lies in using the spatial and temporal variation characteristics of normalized and enhanced vegetation index, not only realizing high-precision inversion of vegetation coverage, but also describing its spatial heterogeneity through probability distribution function, laying a reliable foundation for subsequent analysis. Combined with land surface temperature and thermal inertia model inversion of water content distribution, the vegetation shading effect and soil thermal properties are fully considered, making the water state evaluation more close to the actual conditions. Further, by fusing the coverage time function and water content fluctuation, the dynamic tracking of soil fertility is realized, the first soil fertility time trend reveals the change of fertility from the response of vegetation growth and water stress, and the second soil fertility time trend introduces the physiological time sequence constraint through the abnormal identification of the crop growth period, effectively capturing the growth rhythm abnormality caused by farming activities or environmental stress. Finally, the two types of fertility trends are time-aligned and weighted fused, which not only takes into account the direct coupling effect of vegetation and water, but also integrates the indirect indication of the growth health, so that the generated ecological degradation fluctuation value has both physical mechanism and ecological significance, and can sensitively reflect the gradual degradation or sudden deterioration of cultivated land ecosystem, providing quantitative basis for regional cultivated land resource sustainable management and ecological restoration decision-making.

[0179] The warning module combines the cultivated land area fluctuation value or the ecological degradation fluctuation value of the real cultivated land segmented sub-region with the corresponding preset warning value to perform real-time warning, and automatically generates warning information, and visualizes the change region spatial position, change type and change degree of the real cultivated land segmented sub-region through the GIS platform. It needs to be further explained that the modules in the embodiment include a warning discrimination unit and a real-time labeling unit;

[0180] The warning discrimination unit is based on the cultivated land area fluctuation value and the ecological degradation fluctuation value at consecutive time points. If one of the corresponding cultivated land area fluctuation values or ecological degradation fluctuation values at adjacent two time points does not meet the corresponding warning value, it is determined that there is an abnormality in the real cultivated land segmented sub-region at the corresponding consecutive adjacent time points.

[0181] The real-time labeling unit determines the spatial position and the changed area of the change region according to the absolute value of the time characteristic change value of the area of the adjacent real cultivated land segmentation sub-region determined to exist an abnormality at a continuous time point at the same pixel position point and the first cultivated land area fluctuation value in the corresponding time period and the corresponding pixel position point, determines the type of the corresponding change based on the corresponding normalized vegetation index and enhanced vegetation index in the change region and the corresponding spatial change value, and obtains the change degree of the cultivated land type in the current change region by an evaluation algorithm based on the determined spatial position and the changed area of the change region and the size of the absolute value of the time characteristic change value in the corresponding continuous time period, and displays the spatial position and the changed area, the type of the change and the change degree in real time and dynamically on the real cultivated land segmentation sub-region at the corresponding adjacent time point in the space-time coordinate system through an automatic labeling algorithm.

[0182] It needs to be further explained that the process of abnormality discrimination in the embodiment includes:

[0183] Based on the cultivated land area fluctuation value and the ecological degradation fluctuation value of the real cultivated land segmentation sub-region at the historical continuous time points, the historical maximum value of the cultivated land area fluctuation value and the upper limit of the 95% confidence interval are calculated by a statistical analysis algorithm, the upper limit is determined as the cultivated land area fluctuation warning value, and the area change exceeds the normal range when the value exceeds this value; similarly, the historical maximum value of the ecological degradation fluctuation value and the upper limit of the 95% confidence interval are calculated, the upper limit is determined as the ecological degradation fluctuation warning value, and the ecological degradation degree is abnormal when the value exceeds this value.

[0184] When extracting the fluctuation value of the adjacent time points, based on the cultivated land area fluctuation value sequence and the ecological degradation fluctuation value sequence at the continuous time points, the cultivated land area fluctuation value corresponding to any two adjacent time points is selected by a sequence interception algorithm, denoted as the current area fluctuation value and the previous area fluctuation value, and the corresponding ecological degradation fluctuation value is denoted as the current degradation fluctuation value and the previous degradation fluctuation value, forming a fluctuation value pair of adjacent time points.

[0185] When performing fluctuation value comparison, the current area fluctuation value is compared with the cultivated land area fluctuation warning value by 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. to meet the warning value; 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 respectively to determine whether they meet the warning value.

[0186] When determining the area anomaly, based on the comparison result, a comprehensive determination is made by a logical determination algorithm: if the current area fluctuation value in the two adjacent time points does not satisfy the cultivated land area fluctuation early warning value, that is, is greater than the early warning value, or the previous area fluctuation value does not satisfy the cultivated land area fluctuation early warning value, or the current degradation fluctuation value does not satisfy the ecological degradation fluctuation early warning value, or the previous degradation fluctuation value does not satisfy the ecological degradation fluctuation early warning value, it is determined that the real cultivated land segmentation sub-region under the continuous adjacent time points is abnormal; only when all comparison results satisfy the corresponding early warning value, it is determined to be normal.

[0187] It needs to be further explained that the implementation process of the automatic labeling algorithm in the embodiment includes:

[0188] Based on the adjacent real cultivated land segmentation sub-region determined to exist anomaly under the continuous time points, the absolute value of the time feature change value of each pixel position at the adjacent time points is extracted by a pixel-level change detection algorithm, and the absolute value is compared with a preset fluctuation threshold value to screen out all pixel points exceeding the threshold value; by a spatial clustering analysis algorithm, a DBSCAN clustering method based on density is used, and a neighborhood radius parameter and a minimum pixel number parameter are set, and the pixels with a spatial distance within the neighborhood radius and a pixel number reaching the minimum pixel number requirement are clustered into a cluster to form a continuous change area spatial coordinate set, while the isolated pixel points not meeting the density requirement are removed.

[0189] Based on the change area spatial coordinate set, the actual geographic area corresponding to each pixel is calculated by a pixel area conversion algorithm according to the spatial resolution parameter of the remote sensing image, that is, the actual size of the ground represented by a single pixel.

[0190] The actual geographic areas of all pixels in the change area are accumulated and summed by a region area calculation algorithm to obtain an initial change area; by an area calibration algorithm, an area calibration coefficient is calculated using the first cultivated land area fluctuation value of the adjacent time points, specifically, the first cultivated land area fluctuation value is normalized and added by one as a multiplier, and multiplied by the initial change area to obtain the calibrated actual change area.

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

[0192] By the type determination rule: if the normalized vegetation index difference value and the enhanced vegetation index difference value are both less than zero and the spatial change value is negative, it is determined to be the cultivated land degradation type;

[0193] If the two index difference values are both greater than zero and the spatial change value is positive, it is determined to be the cultivated land recovery type; if the signs of the two index difference values are inconsistent or the absolute value of the spatial change value exceeds the disorder threshold value, it is determined to be the cultivated land utilization transformation type.

[0194] When evaluating the degree of change, the actual change area is divided into three levels of small, medium and large by a hierarchical evaluation algorithm, and the absolute value of the time characteristic change value is divided into three levels of low, medium and high by a hierarchical evaluation algorithm; by a weighted scoring algorithm, the actual change area level is weighted and summed by 60% weight, and the absolute value of the time characteristic change value is weighted and summed by 40% weight, to obtain a comprehensive score value; by a severity division algorithm, the comprehensive score value is divided into three severity levels of slight, moderate and severe according to the three quantile ranges.

[0195] When real-time dynamic labeling is performed, the spatial coordinates of the change area are mapped to the geographic coordinates in the space-time coordinate system by a coordinate conversion algorithm; by a visual labeling algorithm, the cultivated land degradation type area is labeled by a red semi-transparent polygon, the cultivated land recovery type area is labeled by a green semi-transparent polygon, and the cultivated land utilization transformation type area is labeled by a yellow semi-transparent polygon; by a size mapping algorithm, the display size of the labeled polygon is adjusted in proportion to the actual change area value; by a color depth mapping algorithm, the color transparency of the labeled polygon is adjusted according to the severity level, and the transparency of the severe level is the highest, and the transparency of the slight level is the lowest, forming a dynamically updated visual labeling display.

[0196] This process realizes efficient and accurate early warning of abnormal changes in cultivated land resources through the construction of multi-dimensional fluctuation indicators and intelligent discrimination mechanisms, significantly improving the real-time performance and decision support capability of dynamic monitoring of cultivated land. Its core value lies in the comprehensive use of area fluctuation and ecological degradation fluctuation double indicators, which not only captures the sharp changes in the quantity of cultivated land, but also reveals the progressive degradation of the quality of the ecological system, overcoming the one-sidedness of single indicator early warning, making the abnormality determination more comprehensive and reliable. By determining the early warning threshold through historical data statistics, the regional specificity and normal fluctuation range are fully reflected, avoiding false positives and false negatives. Further, the automatic labeling process deeply integrates pixel-level change detection, spatial clustering and multi-source feature analysis, which can accurately locate the spatial range of the change area and quantify the change area, and at the same time, according to the vegetation index change trend and spatial characteristics, intelligently discriminates the change type, clearly distinguishing different situations such as degradation, recovery and utilization transformation. The severity evaluation comprehensively considers the change area and the change intensity, adopts a weighted scoring mechanism, ensuring the scientificity and intuitiveness of the grading results. Finally, through the dynamic visualization labeling on the GIS platform, the abstract data fluctuation is converted into concrete spatial graphics, with color, transparency and size mapping the change type and severity, providing an extremely intuitive and information-rich monitoring view for management personnel.

[0197] Embodiment 2:

[0198] Please refer to Figure 2 Another embodiment provided by the present application: a land ecological condition monitoring method based on remote sensing data, comprising the following steps:

[0199] S1, acquire a remote sensing image sequence of a target area and analyze to obtain an initial cultivated land segmentation sub-region sequence and an edge segmentation feature set;

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

[0201] S3, determine a real cultivated land segmentation sub-region area based on the remote sensing image feature space and the initial cultivated land segmentation sub-region, and determine a cultivated land area fluctuation value based on a determined real cultivated land segmentation sub-region area change value;

[0202] S4, obtain an ecological degradation fluctuation value of the real cultivated land segmentation sub-region based on the remote sensing image feature space, ground surface temperature and the obtained cultivated land vegetation coverage;

[0203] S5, perform real-time early warning based on the cultivated land area fluctuation value or the ecological degradation fluctuation value of the real cultivated land segmentation sub-region and a corresponding preset early warning value, automatically generate early warning information, and visually mark the change area spatial position, change type and change degree of the real cultivated land segmentation sub-region through a GIS platform.

[0204] The embodiments of the present application are described above in combination with the drawings, but the present application is not limited to the above specific embodiments, and the above specific embodiments are only illustrative and not restrictive. Those skilled in the art can make changes, modifications, replacements and variations to the above embodiments without departing from the purpose of the present application and the scope protected by the claims, which are all within the protection of the present application.

Claims

1. A land ecological condition monitoring system based on remote sensing data, characterized in that, The application relates to a remote sensing image-based farmland area and ecological degradation monitoring method and device. The preprocessing module is used for acquiring a remote sensing image sequence of a target area and analyzing and acquiring an initial farmland segmentation sub-region sequence and an edge segmentation feature set; The feature extraction module obtains a remote sensing image feature space based on the initial farmland segmentation sub-region sequence and the edge segmentation feature set; the remote sensing image feature space comprises a spatial feature sequence and a spatial feature fluctuation value of any N edge extension directions of a remote sensing image; the spatial feature sequence of the any N edge extension directions comprises a normalized vegetation index sequence and an enhanced vegetation index sequence corresponding to the extension directions; and the spatial feature fluctuation value comprises a normalized vegetation index spatial change value and an enhanced vegetation index spatial change value of a corresponding remote sensing image at each time point; The region determination module determines a real farmland segmentation sub-region area based on the remote sensing image feature space and the initial farmland segmentation sub-region, and determines a farmland area fluctuation value based on a determined real farmland segmentation sub-region area change value; The inversion module obtains an ecological degradation fluctuation value of the real farmland segmentation sub-region based on the remote sensing image feature space, a ground surface temperature and an inversion obtained farmland vegetation coverage; The warning module performs real-time warning based on the farmland area fluctuation value or the ecological degradation fluctuation value of the real farmland segmentation sub-region and a corresponding preset warning value, simultaneously generates warning information, and visually marks a change region spatial position, a change type and a change degree of the real farmland segmentation sub-region through a GIS platform; The feature extraction module comprises a center determination unit; The center determination unit obtains a first center point based on a pixel position point in the edge segmentation feature set corresponding to the initial farmland segmentation sub-region at each time point and a minimum circumscribed rectangle, obtains a centroid position as a second center point based on the edge pixel position point, and determines a real center position point of the initial farmland segmentation sub-region at each time point based on a mean value of the first center point and the second center point; A two-dimensional coordinate is constructed based on the real center position point at each time point in the initial farmland segmentation sub-region at a corresponding time, a line of all aligned real center position points at continuous time points is taken as a three-dimensional time axis, a space-time coordinate system is constructed, and the initial farmland segmentation sub-region of the remote sensing image collected at each time point is embedded into the space-time coordinate system at a corresponding time point; The region determination module comprises a region area unit; The region area unit obtains a real farmland segmentation sub-region area at each time point based on a second edge position point set of the initial farmland segmentation sub-region at each time point and the corresponding initial farmland segmentation sub-region and a region integration algorithm, and takes each position pixel point of the real farmland segmentation sub-region with the largest area as a starting point to extract an absolute value of a time feature change value of a same pixel position point of adjacent real farmland segmentation sub-region areas at continuous time points along positive and negative axes of the three-dimensional time axis. The preprocessing module comprises an acquisition unit and an image segmentation unit; 2. The land ecological condition monitoring system based on remote sensing data according to claim 1, characterized in that, ​ The acquisition unit is configured to acquire a remote sensing image sequence of a target region in a continuous preset time length, and perform radiation correction and geometric correction preprocessing on the remote sensing image sequence to obtain a corrected remote sensing image sequence. The image segmentation unit is configured to perform initial region identification on the corrected remote sensing image sequence in combination with a pre-trained YOLOv5 model, to obtain an initial region identification range and a corresponding region type of the remote sensing image collected at different time points, and to obtain an edge segmentation feature set corresponding to each initial farmland segmentation sub-region at each time point based on the initial region identification range and an image segmentation algorithm. The region type includes farmland, construction land, woodland, pond, and abandoned land.

3. The land ecological condition monitoring system based on remote sensing data according to claim 2, characterized in that, The feature extraction module further includes a feature extraction unit. The feature extraction unit is configured to take a real center position point of a corresponding initial farmland segmentation sub-region at each time point as a starting point, and perform feature extraction along any N edge extension directions of the initial farmland segmentation sub-region at the corresponding time point, to obtain a spatial feature sequence and a spatial feature fluctuation value of any N edge extension directions at the current time point.

4. The land ecological condition monitoring system based on remote sensing data according to claim 3, characterized in that, The region determination module further includes a region determination unit. The region determination unit is configured to combine a region type fluctuation threshold with the spatial feature sequence and the spatial feature fluctuation value of any N edge extension directions at the current time point, and when the spatial feature fluctuation value of M consecutive pixel position points of any edge extension direction does not satisfy the region type fluctuation threshold, to take a pixel position point corresponding to a maximum spatial feature fluctuation value as an edge point of the current edge extension direction, and repeat the operation N times to determine a second edge position point set.

5. The land ecological condition monitoring system based on remote sensing data according to claim 4, characterized in that, The region determination module further includes a region fluctuation determination unit. The region fluctuation determination unit is configured to obtain a first farmland area fluctuation value sequence according to areas of real farmland segmentation sub-regions corresponding to adjacent time points, and obtain a proportion of absolute values of time feature change values of the same pixel position points that do not satisfy a region type fluctuation threshold as a second farmland area fluctuation value sequence based on the absolute values of the time feature change values and the region type fluctuation threshold. The region fluctuation determination unit is further configured to obtain a farmland area fluctuation value of an adjacent time point based on a mean value of a first farmland area fluctuation value and a second farmland area fluctuation value of the corresponding adjacent time points in the first farmland area fluctuation value sequence and the second farmland area fluctuation value sequence.

6. The land ecological condition monitoring system based on remote sensing data according to claim 5, wherein, The inversion module includes a coverage inversion unit. The coverage inversion unit is configured to obtain a farmland vegetation coverage and a coverage distribution probability function corresponding to each real farmland segmentation sub-region at a corresponding time point based on a normalized vegetation index sequence and an enhanced vegetation index sequence corresponding to each time point, and the normalized vegetation index spatial change value and the enhanced vegetation index spatial change value at the corresponding time point, and to obtain a farmland vegetation coverage time distribution function based on the farmland vegetation coverage corresponding to each real farmland segmentation sub-region and the absolute values of time feature change values of adjacent real farmland segmentation sub-region areas at the same pixel position point at consecutive time points.

7. The land ecological condition monitoring system based on remote sensing data according to claim 6, characterized in that, The inversion module further includes a first change trend unit. The first change trend unit is configured for: based on the corresponding farmland vegetation coverage and the corresponding coverage distribution probability function of each real farmland segmentation sub-region, combining the ground surface temperature and the thermal inertia model, inversely deriving the water content distribution state function of the corresponding real farmland segmentation sub-region; and based on the water content distribution state function of the real farmland segmentation sub-region, combining the farmland vegetation coverage time distribution function, obtaining the water content time fluctuation function of the corresponding real farmland segmentation sub-region at continuous time points. Based on the corresponding farmland vegetation coverage and the corresponding coverage distribution probability function of each real farmland segmentation sub-region, combining the water content distribution state function of each real farmland segmentation sub-region, the soil fertility distribution function of each real farmland segmentation sub-region is inversely derived; and based on the soil fertility distribution function of each real farmland segmentation sub-region at continuous time points and the farmland vegetation coverage time distribution function, the first soil fertility time change trend at continuous time points is inversely derived.

8. The land ecological condition monitoring system based on remote sensing data according to claim 7, characterized in that, The inversion module further comprises a second change trend unit and a comprehensive trend unit. The second change trend unit is configured for: based on the normalized vegetation index spatial change value and the enhanced vegetation index spatial change value, combining the crop phenology extraction algorithm, obtaining the corresponding crop phenology node of the real farmland segmentation sub-region; based on the time characteristic change value corresponding to the corresponding crop phenology node of the real farmland segmentation sub-region, performing phenology anomaly identification, obtaining the corresponding crop phenology anomaly identification characteristic space; and based on the corresponding crop phenology anomaly identification characteristic space, inversely deriving the second soil fertility time change trend of the corresponding real farmland segmentation sub-region. The comprehensive trend unit is configured for: performing time alignment processing on the first soil fertility time change trend and the second soil fertility time change trend, and using the first soil fertility time change trend and the second soil fertility time change trend in the aligned time period to obtain the ecological degradation fluctuation value of the real farmland segmentation sub-region by combining a weighted average algorithm.

9. The land ecological condition monitoring system based on remote sensing data according to claim 8, wherein, The early warning module comprises a warning discrimination unit and a real-time labeling unit. The warning discrimination unit is configured for: based on the farmland area fluctuation value and the ecological degradation fluctuation value at continuous time points, if one of the corresponding farmland area fluctuation values or the ecological degradation fluctuation values at adjacent two time points does not satisfy the corresponding early warning value, it is determined that there is an anomaly in the real farmland segmentation sub-region at the corresponding continuous adjacent time points. The real-time labeling unit determines the spatial position and the changed area of the change region according to the absolute value of the time feature change value of the adjacent real cultivated land segmentation sub-region area at the continuous time points at the same pixel position point and the first cultivated land area fluctuation value and the corresponding pixel position point in the corresponding time period, determines the corresponding change type based on the corresponding normalized vegetation index and enhanced vegetation index in the change region and the corresponding spatial change value, and obtains the change degree of the cultivated land type in the current change region through an evaluation algorithm based on the determined spatial position and the changed area of the change region and the size of the absolute value of the time feature change value in the corresponding continuous time period, and displays the spatial position and the changed area, the change type and the change degree in real time and dynamically on the real cultivated land segmentation sub-region corresponding to the adjacent time points in the spatio-temporal coordinate system through an automatic labeling algorithm.

10. A method for monitoring land ecological conditions based on remote sensing data, which is implemented based on the system for monitoring land ecological conditions based on remote sensing data according to any one of claims 1-9, characterized in that, It comprises: acquiring a remote sensing image sequence of a target region and analyzing to obtain an initial cultivated land segmentation sub-region sequence and an edge segmentation feature set; based on the initial cultivated land segmentation sub-region sequence and the edge segmentation feature set, obtaining a remote sensing image feature space; based on the remote sensing image feature space and the initial cultivated land segmentation sub-region, determining a real cultivated land segmentation sub-region area, and determining a cultivated land area fluctuation value based on the determined real cultivated land segmentation sub-region area change value; based on the remote sensing image feature space, the ground surface temperature and the obtained cultivated land vegetation coverage, obtaining an ecological degradation fluctuation value of the real cultivated land segmentation sub-region; based on the cultivated land area fluctuation value or the ecological degradation fluctuation value of the real cultivated land segmentation sub-region and the corresponding preset warning value, performing real-time warning, automatically generating warning information, and visualizing and labeling the change region spatial position, change type and change degree of the real cultivated land segmentation sub-region through a GIS platform.

Images