A method and system for remote sensing monitoring of cultivated land protection

Abstract

The present application relates to the technical field of image processing, in particular to a remote sensing monitoring method and system for cultivated land protection, comprising the following steps: acquiring a remote sensing image sequence of a cultivated land area at different time nodes, extracting the pixel gray value of the same land block and arranging it according to time, calculating the ratio of the difference value between adjacent frames and the gray mean value in the frame, marking the position continuously exceeding the fluctuation reference ratio, and generating a pixel region set. In the present application, by adopting a sampling frequency threshold optimization strategy, combining the albedo variation rate to screen abnormal areas, reducing short-term fluctuation interference, enhancing time sequence consistency, introducing the solar incident angle and terrain slope calculation projection influence to eliminate the slope shielding area, using spatial grid marking and boundary integrity judgment to improve the coherence of cultivated land distribution information, combining time sequence analysis, sampling optimization, physical property compensation and spatial consistency correction to enhance the fine identification of cultivated land dynamic change, and improving the monitoring accuracy and data reliability.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of image processing, and in particular to a remote sensing monitoring method and system for cultivated land protection. BACKGROUND

[0002] The technical field of image processing includes related technologies for collecting, analyzing, recognizing, converting, and storing images using computers or other electronic devices. The core content of this field includes image acquisition, image preprocessing, feature extraction, image classification and recognition, image compression and coding, etc. Image processing is widely used in medical image analysis, industrial detection, remote sensing monitoring, security monitoring, autonomous driving, and other fields, and its overall development relies on the continuous improvement of image sensor technology, image computing models, and image data processing methods. This field focuses on improving image clarity and enhancing the amount of useful information in images, and achieving effective analysis and interpretation of image content through mathematical modeling and computing technology.

[0003] Among them, the remote sensing monitoring method for cultivated land protection refers to a method of observing and analyzing the range and changes of cultivated land through remote sensing images. This patent subject addresses the problem of remote sensing identification in cultivated land protection, covering techniques such as obtaining cultivated land distribution information through multi-temporal remote sensing images, identifying cultivated land change areas using image difference analysis methods, and determining cultivated land properties in combination with vegetation indices and land use types. This method generally uses image contrast analysis, surface coverage classification determination, and time series image change identification to complete remote sensing monitoring operations, thereby achieving accurate extraction and change monitoring of cultivated land information at the data level.

[0004] The existing technology relies on fixed time node image comparison in cultivated land remote sensing monitoring, and does not fully utilize pixel-level gray sequence information, leading to misjudgment caused by short-term environmental disturbance, affecting the stable identification of change areas. The sampling method is not adjusted according to the change frequency, resulting in data redundancy in some areas and information loss in some areas, affecting the balance of overall monitoring. Time series monitoring fails to effectively eliminate short-term abnormal fluctuations, leading to a decrease in local data reliability and affecting the accurate judgment of cultivated land change trends. Slope shading is not effectively identified, and some areas lack surface information due to shadow or reflection errors, affecting the accurate delineation of cultivated land boundaries. The spatial information expression method does not fully consider the integrity of plot boundaries, resulting in broken information in some cultivated land plots, which is not conducive to the coherent analysis of large-scale cultivated land changes. These shortcomings affect data precision and consistency, reducing the reliability of cultivated land monitoring results. SUMMARY

[0005] The purpose of the present application is to solve the shortcomings in the prior art and to provide a remote sensing monitoring method and system for cultivated land protection.

[0006] In order to achieve the above object, the present application adopts the following technical scheme: A remote sensing monitoring method for cultivated land protection, comprising the following steps:

[0007] S1: acquiring a remote sensing image sequence of the cultivated land region at different time nodes, extracting pixel gray values of the same land block and arranging them in time, calculating the ratio of the difference value between adjacent frames and the gray mean value in the frame, marking the positions continuously exceeding the fluctuation reference ratio, and generating a pixel region set;

[0008] S2: calling the pixel region set, setting an interference frequency threshold to screen frequent regions, and checking the position and coverage relationship frame by frame, performing priority reduction on the pixels in the sampling dense region, and generating a time sequence image sampling priority reference value;

[0009] S3: calling the time sequence image sampling priority reference value, extracting a multi-band albedo sequence of the corresponding region, setting a time window and calculating the inter-frame variation rate, comparing the current frame variation rate with the fluctuation value in the window, and obtaining a sampling reserved region segment set;

[0010] S4: calling the sampling reserved region segment set, extracting the solar incident angle and the slope parameter, calculating the illumination projection length and combining the reflection intensity change amplitude, screening the regions with combined values exceeding the shielding threshold, and obtaining a shielding excluded image position set;

[0011] S5: calling the shielding excluded image position set, marking the remaining region, tracking the extension state of the land block boundary, judging the boundary integrity and continuity, and generating a sampleable cultivated land image block distribution information.

[0012] As a further scheme of the present application, the pixel region set includes pixel position, gray fluctuation ratio, change continuity mark, the time sequence image sampling priority reference value includes sampling region, sampling frequency threshold, pixel sampling order, the sampling reserved region segment set includes albedo sequence, fluctuation rate, jump threshold, the shielding excluded image position set includes solar incident angle, terrain slope, projection length, shielding determination threshold, and the sampleable cultivated land image block distribution information includes spatial index, position set, and sampling continuous region.

[0013] As a further scheme of the present application, the specific steps of S1 are as follows:

[0014] S101: acquiring a remote sensing image sequence of the cultivated land region at different time nodes, extracting pixel gray values of the same region and sorting them, calculating the gray difference value of adjacent time nodes, combining the gray mean value in the image frame, calculating the ratio of the gray difference value of the pixel points, and obtaining a pixel gray ratio sequence;

[0015] S102: Screen the pixel positions whose gray scale ratio changes more than the set reference ratio at continuous time nodes based on the pixel gray scale ratio sequence, and collect the gray scale fluctuation marked pixel point set after screening;

[0016] S103: Call the gray scale fluctuation marked pixel point set, analyze the spatial distribution characteristics of adjacent pixels, merge high fluctuation pixel points, and form a pixel region set.

[0017] As a further scheme of the present application, the gray scale difference value calculation formula of the adjacent time nodes is specifically:

[0018] ;

[0019] Among them, represents the gray scale difference value of the pixel point (x, y) between the i-th and i+1-th frames of images, represents the gray scale value of the pixel point (x, y) in the i-th frame, represents the average value of the gray scale values of all pixel points in the i-th frame image, represents the total number of pixels in the image, represents the algebraic sum of the gray scale value changes of all pixel points in the image between the i-th and i+1-th frames.

[0020] As a further scheme of the present application, the specific steps of S2 are:

[0021] S201: Call the pixel region set, set a pixel frequency threshold, mark the regions exceeding the threshold, check the spatial position of the marked regions frame by frame, analyze the coverage range, set the threshold and select the regions with higher frequency, and obtain the high-frequency marked region distribution;

[0022] S202: Based on the high-frequency marked region distribution, perform sampling and reduction operation at the selected position in the sampling dense region, while retaining the original sampling data level of the unmarked region, and obtain the adjusted pixel sampling level;

[0023] S203: Call the adjusted pixel sampling level, record the sampling order structure of the pixel points, arrange in time sequence, and establish the time sequence image sampling priority reference value.

[0024] As a further scheme of the present application, the specific steps of S3 are:

[0025] S301: Call the time sequence image sampling priority reference value, extract the multi-band albedo sequence corresponding to the adjusted sampling region, set a time window and calculate the albedo fluctuation rate between continuous frames, arrange the fluctuation rate of the time nodes, and establish the albedo fluctuation rate sequence;

[0026] S302: Based on the albedo variation rate sequence, a difference between the variation rate of the region in the current sampling frame and the fluctuation rate in the time window is compared, regions with a difference exceeding an albedo jump threshold are excluded, and an albedo stable region is obtained;

[0027] S303: The albedo stable region is called to extract spatial distribution information of the reserved region, region segment data is arranged in a time sequence, and a sampling reserved region segment set is generated.

[0028] As a further scheme of the present application, the albedo variation rate calculation formula is specifically:

[0029] ;

[0030] Among them, represents the albedo variation rate of the pixel point (x, y) in the i-th image in the wave band b, represents the albedo value of the pixel point (x, y) in the i-th image in the wave band b, represents the average albedo value of all pixel points in the wave band b in the frame image, represents the albedo value of the pixel point (x, y) in the selected reference time frame in the wave band b, represents the time interval between the current image frame and the reference image frame.

[0031] As a further scheme of the present application, the specific steps of S4 are:

[0032] S401: The sampling reserved region segment set is called to extract the solar incident angle and the terrain slope value of the corresponding region, the projection length in the solar direction is calculated based on the incident angle and the slope value, the illumination projection range of the region is identified, and whether there is a slope shielding condition is judged, and a slope shielding judgment result is obtained;

[0033] S402: Based on the slope shielding judgment result, the reflection intensity variation amplitude of the region is extracted, and combined operation is performed with the projection value, the combined value of the region is calculated, and a shielding influence combined value is obtained;

[0034] S403: The shielding influence combined value is called to screen regions with a combined value exceeding a shielding judgment threshold, extract corresponding image positions and arrange spatial information, and generate a shielding excluded image position set.

[0035] As a further scheme of the present application, the specific steps of S5 are:

[0036] S501: Call the shielding exclusion image position set, mark the distribution state of the remaining area in the spatial grid, track the extension form of the corresponding land boundary in the image frame, analyze the spatial connection of the area according to the continuity characteristics of the boundary adjacent area, and judge the integrity and continuity of the boundary to obtain a boundary coherent area;

[0037] S502: Based on the boundary coherent area, filter the area meeting the sampling structure continuity judgment condition, classify the area meeting the condition as a candidate block, and obtain a sampleable candidate area;

[0038] S503: Call the sampleable candidate area, organize the spatial index information, and match the corresponding image position to establish the sampleable cultivated land image block distribution information.

[0039] A remote sensing monitoring system for cultivated land protection comprises:

[0040] The time sequence gray scale change extraction module obtains a remote sensing image sequence of the cultivated land area at different time nodes, extracts the gray scale values of the same block of land pixels and sorts them, calculates the gray scale difference value between adjacent time nodes, judges the difference value by combining the gray scale mean value in the image frame, marks the pixels exceeding the gray scale fluctuation reference value in continuous multiple changes, and generates a pixel area set;

[0041] The high-frequency variable region screening module calls the pixel area set, sets a frequency threshold and labels a high-frequency region, checks the spatial position and coverage relationship of the labeled region, adjusts the sampling level of the dense sampling region, records the adjusted pixel sampling order, and generates a time sequence image sampling priority reference value;

[0042] The albedo fluctuation screening module calls the time sequence image sampling priority reference value, extracts the multi-band albedo sequence of the adjusted sampling area, calculates the albedo change rate in the time window, screens out the regions whose rate difference exceeds the albedo jump threshold, and obtains a sample retention region segment set;

[0043] The slope shielding exclusion module calls the sample retention region segment set, extracts the solar incident angle and terrain slope value of the corresponding region, calculates the solar direction projection length, judges the shielding condition, combines the reflection intensity variation amplitude and the projection value of the region, screens out the regions exceeding the shielding judgment threshold, and obtains a shielding exclusion image position set;

[0044] The cultivated land image block extraction module calls the shielding exclusion image position set, marks the spatial distribution of the remaining area, tracks the land boundary in the image frame, judges the integrity and continuity of the boundary, screens the area meeting the sampling structure continuity, outputs the spatial index and position set, and generates the sampleable cultivated land image block distribution information.

[0045] Compared with the prior art, the advantages and positive effects of the present application are that:

[0046] In the present application, by adopting the sampling frequency threshold optimization strategy, combining the albedo variation rate screening abnormal area, reducing short-term fluctuation interference, enhancing the time sequence consistency, introducing the solar incident angle and terrain slope calculation projection influence, eliminating the slope shielding area, adopting the spatial grid marking and boundary integrity judgment, improving the continuity of the cultivated land distribution information, combining the time sequence analysis, sampling optimization, physical property compensation and spatial consistency correction, enhancing the fine identification of the cultivated land dynamic change, improving the monitoring precision and data reliability. BRIEF DESCRIPTION OF DRAWINGS

[0047] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed to be used in the embodiment description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0048] Figure 1 The step flowchart of the present application is shown in the figure.

[0049] Figure 2 The system module diagram of the present application is shown in the figure. DETAILED DESCRIPTION

[0050] The technical solutions in the present application will be described below in combination with the drawings.

[0051] In the embodiments of the present application, the words such as "example", "for example" are used to represent as an example, illustration or description. Any embodiment or design scheme described as "example" in the present application should not be interpreted as more preferred or more advantageous than other embodiments or design schemes. Rather, the word "example" is intended to present the concept in a specific way. In addition, in the embodiments of the present application, the meaning expressed by "and / or" can be both, or can be one of the two.

[0052] In the embodiments of the present application, "image" and "picture" can be used interchangeably at times. It should be pointed out that when the distinction is not emphasized, the meanings expressed are consistent. "Of", "corresponding" and "corresponding" can be used interchangeably at times. It should be pointed out that when the distinction is not emphasized, the meanings expressed are consistent.

[0053] In the embodiments of the present application, sometimes the subscript such as W1 can be written in the form of non-subscript such as W1. When the distinction is not emphasized, the meanings expressed are consistent.

[0054] In order to make the technical problems, technical solutions and advantages to be solved by the present application more clear, the following will be described in detail in combination with the drawings and specific embodiments.

[0055] Please refer to Figure 1 A remote sensing monitoring method for cultivated land protection, comprising the following steps:

[0056] S1: Obtain a remote sensing image sequence of the cultivated land area at different time nodes, extract the gray values of the same block of pixels in the image and arrange them in time sequence, perform point-by-point operation on the gray value difference between adjacent time nodes, combine the gray value difference with the gray mean value in the image frame to perform ratio judgment, mark the pixel positions that exceed the gray fluctuation reference ratio in continuous multiple changes, and generate a pixel region set;

[0057] S2: Call the pixel region set, set a frequency threshold and mark the regions with high frequency, check the spatial position and coverage relationship of the marked regions in the image frame by frame, perform sampling and reduction operation on the positions in the dense region, while retaining the original sampling level of the unmarked region, record the adjusted pixel sampling order structure, and generate a time sequence image sampling priority reference value;

[0058] S3: Call the time sequence image sampling priority reference value, extract the multi-band albedo sequence corresponding to the adjusted sampling region in the image frame, define a time window and calculate the albedo variation rate between consecutive frames in the sequence, compare the variation rate of each region in the current sampling frame with the difference value of the fluctuation rate in the time window, exclude the regions whose difference value exceeds the albedo jump threshold, and obtain a sampling reserved region segment set;

[0059] S4: Call the sampling reserved region segment set, extract the solar incident angle and terrain slope value of the corresponding region in the image, calculate the projection length in the solar direction for each region, judge whether there is a slope shielding condition, and then select the regions whose combined value exceeds the shielding determination threshold according to the variation amplitude of the region reflection intensity and the projection value, to obtain a shielding excluded image position set;

[0060] S5: Call the shielding excluded image position set, mark the distribution state of the remaining regions in the spatial grid, track the corresponding land boundary extension form in the image frame, perform boundary integrity and continuity judgment on the regions, select the regions that meet the sampling structure continuous determination condition as candidate blocks, and output the spatial index and position set, to generate the sampleable cultivated land image block distribution information.

[0061] The pixel region set comprises a pixel position, a gray scale fluctuation ratio, and a change continuity mark. The time sequence image sampling priority reference value comprises a sampling region, a sampling frequency threshold, and a pixel sampling sequence. The sampling reserved region segment set comprises an albedo sequence, a fluctuation rate, and a jump threshold. The shielding exclusion image position set comprises a solar incident angle, a terrain slope, a projection length, and a shielding determination threshold. The cultivable land image block distribution information comprises a spatial index, a position set, and a sampling coherent region.

[0062] The specific steps of S1 are as follows:

[0063] S101: Obtain a remote sensing image sequence of a cultivable land region at different time nodes, extract pixel gray scale values of the same region and sort them, calculate a gray scale difference value between adjacent time nodes, combine an image frame gray scale mean value, perform ratio calculation on the pixel gray scale difference value, and obtain a pixel gray scale ratio sequence.

[0064] The gray scale difference value calculation formula between adjacent time nodes is as follows:

[0065] ;

[0066] Wherein, represents a gray scale difference value of a pixel point (x, y) between the i-th and the i+1-th frames of images, represents a gray scale value of the pixel point (x, y) in the i-th frame, represents an average value of all pixel gray scale values in the i-th frame of image, represents a total number of pixels in the image, represents an algebraic sum of gray scale value changes of all pixel points in the image between the i-th and the i+1-th frames;

[0067] The formula is used to calculate the gray scale difference value of a pixel point in a remote sensing image between adjacent time nodes, and integrates local gray scale difference, image overall gray scale change trend and balanced correction term of all pixel changes. The parameters are explained as follows:

[0068] The image frames T1 and T2 correspond to remote sensing images obtained on April 5, 2024 and May 5, 2024 respectively. The image size is 200 rows x 200 columns, a total of 40000 pixel points. The image type is a single-channel gray scale image, and the gray scale value range is 0 to 255. The gray scale values of the pixel point (x=100, y=100) in the T1 frame and the T2 frame are obtained by image reading, and the gray scale values are 122 and 139 respectively, that is:

[0069] ;

[0070] The average value of all pixel gray scale values in the image frame T1 is obtained by summing all pixels and then dividing by the total number of pixels. The total sum is 4928000, and the total number of image pixels is 40000 pixel points:

[0071] ;

[0072] The total sum of pixel gray scale in image frame T2 is 5260000, and the calculation is as follows:

[0073] ;

[0074] The sum of the gray scale changes of all pixels in image frames T1 to T2 is the total difference value monitored by the system, which is 212000, and the image pixel number N = 40000 is substituted:

[0075] ;

[0076] The above items are substituted into the formula:

[0077] The first term is the square difference of the gray scale change of the pixel point itself:

[0078] ;

[0079] The second term is the square difference of the mean value change of the image frame:

[0080] ;

[0081] After adding the two terms, the root value is obtained:

[0082] ;

[0083] The third term is the absolute value of the average difference of the pixel full image change:

[0084] ;

[0085] The final result is:

[0086] ;

[0087] The result shows that the pixel point (100, 100) has a gray scale difference of 13.61 compared to the overall image change trend between April and May 2024, representing a moderate brightening trend of the pixel point in local and global gray scale changes. This value is used in subsequent calculations of image gray scale mean value ratio to construct a pixel gray scale ratio sequence. The larger the value, the more unstable the point is in the local interframe change, and the more significant the difference with the overall image.

[0088] S102: Based on the pixel gray scale ratio sequence, filter the pixel positions whose gray scale ratio changes at consecutive time nodes exceed a set reference ratio, and collect the filtered pixel positions to obtain a set of gray scale fluctuation marked pixel points;

[0089] The ratio of each pixel at consecutive time nodes is compared, and the pixel points with obvious gray scale changes at two or more consecutive time nodes are screened out. To define "obvious gray scale change", a ratio threshold value needs to be set, for example, referring to the median and standard deviation of each ratio under large-scale statistics, the threshold value is set to 0.15, and all pixel ratio sequences are traversed. For each pixel, it is judged whether the ratio of any two consecutive nodes is greater than 0.15, for example, the ratio sequence of pixel (100, 100) is 0.117, 0.167, 0.138. Since the second ratio is greater than 0.15 and the third ratio is slightly less than 0.15, although there is a node that exceeds the threshold value, it does not meet the condition of more than two consecutive nodes greater than the threshold value, so it is not selected. Looking at the ratio sequence of pixel (120, 120), which is 0.162, 0.185, and 0.190, all of which are greater than 0.15, and are consecutive time nodes, which meet the conditions, the pixel is marked as a gray scale fluctuation pixel and stored in the marked set. Using a two-dimensional coordinate array form, assuming that the total number of image pixels is 40,000, about 1,200 pixels meet the conditions, forming a gray scale fluctuation marked pixel point set. In practical applications, such fluctuation pixel points often concentrate in the image nodes corresponding to the rapid growth or harvesting period of crops, corresponding to the crop replacement or management behavior in the farmland, and have actual reference to the entire image scene.

[0090] S103: Call the gray scale fluctuation marked pixel point set, analyze the spatial distribution characteristics of adjacent pixels, merge high fluctuation pixel points, and form a pixel region set.

[0091] Further identification of its spatial distribution is needed to determine whether a contiguous area is formed. The method is to traverse each pixel position in the marked set one by one, check the positions of the eight adjacent pixel points around it, and confirm whether they also exist in the gray scale fluctuation set. For example, for pixel (120, 120), check whether at least three of the eight positions (119, 119), (119, 120), (119, 121), (120, 119), (120, 121), (121, 119), (121, 120), (121, 121) around it are also in the set. If the condition is met, the pixel and its neighbors are considered as a preliminary aggregation unit, and multiple pixels that meet the condition are aggregated and grouped into the same group to form a larger area. Through horizontal and vertical scanning of the image, these adjacent relationships are continuously searched, and each connected region is assigned a unique number to form a pixel region set. For example, six change areas are finally obtained, each containing 85, 74, 60, 52, 39 and 25 pixels, respectively, accounting for 0.42%, 0.37%, 0.3%, 0.26%, 0.2% and 0.13% of the image area, respectively. The first two regions show a rectangular aggregation trend, which may reflect the crop replacement process within a large cultivated land, and the last few regions show a scattered dot aggregation, which may be related to local management of farmland. The entire region set finally forms a pixel region set, which serves as the basis data for subsequent analysis of regional change trends.

[0092] The specific steps of S2 are:

[0093] S201: Call the pixel region set, set the pixel frequency threshold, mark the regions exceeding the threshold, check the spatial position of the marked regions frame by frame, analyze the coverage, set the threshold and select the regions with higher frequency, and obtain the distribution of high-frequency marked regions.

[0094] In the calling stage, the set of regions with completed grayscale change annotation is processed, and the number of times each region is marked as a changed region in all time nodes is calculated, that is, the frequency of each region being identified as a changed region in the image frame sequence is counted. For example, region Z1 appears 4 times in 5 images, and its frequency is 4. When setting the pixel frequency threshold, the overall distribution of all region frequencies should be considered. The maximum value, minimum value and average value of all region frequencies can be counted first. For example, all region frequency values are between 1 and 6, and the average value is 3.2. The frequency threshold can be set to an integer greater than the average value, that is, set to 4, which means that only when a region is labeled as a changed region at least 4 times in an image sequence, the region is considered to have sufficient stable change characteristics. Regions with a frequency greater than or equal to 4 are labeled. The pixel position of the labeled region in the image is recorded as a two-dimensional coordinate, for example, region Z1 contains pixel positions (120, 120), (121, 120), (120, 121), and (121, 121). In each frame, check whether these positions are within the image boundary range and have annotation records. If there are offset or discontinuous regions, mark them as position drift regions. Further, calculate the pixel coverage of each high-frequency region in each frame, that is, calculate the number of all labeled pixels in the region, which is used as an area statistical indicator. Compare the coverage ranges of all high-frequency regions, and select regions with an area greater than a set lower limit, for example, set the minimum pixel coverage to 50 pixels. If a region only covers 40 pixels, it is excluded. Finally, the regions that meet the conditions in frequency and area are high-frequency labeled regions. The spatial distribution of these regions is then counted and sorted, for example, regions Z1, Z3 and Z5 are classified into the same position band to form a group of high-frequency distribution regions. Finally, the high-frequency labeled region distribution is output.

[0095] S202: Based on the high-frequency labeled region distribution, perform sampling and downscaling operations in the selected position of the sampling dense region, while retaining the original sampling data level of the unlabeled region, and obtain the adjusted pixel sampling level.

[0096] First, it needs to be confirmed which areas are in the image sampling dense area, the identification of sampling dense area is completed by statistics of the number of sampling points per unit area in each frame of the image, for example, the image is divided into 10x10 grid blocks, the number of labeled regions and sampling points in each grid block is counted, if the coverage ratio of high-frequency labeled regions in a certain grid is greater than 80%, and the average number of pixels in each block contains more than 20 sampling points, then the grid block is determined as a sampling dense area, record the spatial position and corresponding pixel index of all such areas, perform sampling reduction on these positions, that is, reduce the sampling points per unit area to a certain proportion of the original, for example, adjust from one sampling point per 64 pixels to one sampling point per 128 pixels, that is, originally a certain area in a frame of image contains 20 sampling points, after adjustment, only 10 points are retained, and the reserved strategy preferentially retains the boundary pixels and the center pixels of the region, and deletes the sampling points in the middle overlapping area. The sampling reduction operation is recorded in the pixel sampling level matrix. For non-high-frequency labeled regions, that is, unlabeled regions, do not change the sampling level, keep the original sampling structure, for example, a certain region has one sampling point per 128 pixels in the original image, keep this ratio unchanged. Finally, integrate the sampling level data of different regions to obtain the complete adjusted pixel sampling level matrix containing high-frequency regions and unlabeled regions.

[0097] S203: Call the adjusted pixel sampling level, record the sampling order structure of the pixel points, arrange in time sequence, and establish the time sequence image sampling priority reference value;

[0098] The structure of the matrix needs to be arranged, that is, the sampling level value corresponding to each pixel point is sorted and arranged in the time sequence of the image frame, for example, the sampling level value of pixel (100, 100) in image frames T1 to T6 is 1, 1, 2, 1, 1, 2, respectively, indicating that it is downgraded in the 3rd and 6th frames. According to the time sequence, the sampling change trajectory of the pixel is constructed, and the sampling level of all pixel points is arranged in the time axis sequence to generate a sampling order structure table, each row is a pixel point, and each column is the sampling level value of the corresponding frame. In the process of establishing the time sequence image sampling priority, the sampling level sequence of each pixel point is analyzed, if it remains low level sampling in most time frames, it is marked as a low priority pixel point, if it remains high level sampling in multiple time frames, it is marked as a high priority, for example, pixel point (150, 150) remains high sampling level in 5 frames out of 6 frames, then it is arranged in the high priority sequence, and finally a sampling priority matrix arranged in time frames is formed, which records the sampling change trend of different pixels in the image sequence. The constructed sampling priority structure serves as the basis for subsequent image data reading or processing.

[0099] The specific steps of S3 are:

[0100] S301: Call the time sequence image sampling priority reference value, extract the multi-band albedo sequence corresponding to the adjusted sampling area, set the time window and calculate the albedo variation rate, arrange the variation rate for the time node, and establish the albedo variation rate sequence;

[0101] The albedo variation rate calculation formula is specifically:

[0102] ;

[0103] Wherein, represents the albedo variation rate of pixel point (x, y) in the i-th frame image in band b, represents the albedo value of pixel point (x, y) in the i-th frame image in band b, represents the average albedo value of all pixel points in the frame image in band b, represents the albedo value of pixel point (x, y) in the selected reference time frame in band b, represents the time interval between the current image frame and the reference image frame;

[0104] The formula is used to quantify the local variation rate of the albedo of the pixel point in the remote sensing image in a time frame i relative to the reference frame. Select the land parcel number HJ-A12, and the center pixel coordinates of the corresponding area are (160, 160). The albedo information of the pixel in the red band is extracted. According to the Sentinel-2 image data obtained on June 5, 2024, the albedo value after atmospheric correction is read as

[0105] ;

[0106] The image mean value of the red band in the image is obtained by dividing the sum of all effective pixel albedo values by the number of effective pixels, the total albedo value is 201835, and the number of pixels is 62500 pixels, and the calculation result is:

[0107] ;

[0108] The reference time frame is set as the image obtained on May 20, 2024, which is an image after sowing, and is used as the initial reference frame for judging the change of ground albedo. The albedo value of the pixel in the red band is:

[0109] ;

[0110] The image shooting date interval is obtained by subtracting the time stamp in the metadata, and the interval between June 5, 2024 and May 20, 2024 is 16 days, and the result is:

[0111] ;

[0112] Substitute all values into the formula to calculate:

[0113] The first term is the deviation of the current frame from the image mean:

[0114] ;

[0115] The second term is the albedo difference between the current frame and the reference frame:

[0116] ;

[0117] The numerator is the sum of the two:

[0118] ;

[0119] The final rate of change is:

[0120] ;

[0121] This result shows that the albedo rate of change of pixel point (160, 160) in the red band is 0.0024375, with a unit of albedo value divided by days, indicating that within a 16-day interval, the albedo of this pixel changes by about 0.0024 per day on average. This value is used to screen stable or abnormal ground object areas when compared with other pixels or regions in subsequent comparisons. By sequentially arranging with other bands and pixels, a complete albedo rate of change sequence can be formed for time series remote sensing analysis.

[0122] S302: Based on the albedo rate of change sequence, compare the difference between the rate of change of the region in the current sampling frame and the rate of change within the time window, and screen out regions with a difference exceeding the albedo jump threshold to obtain stable albedo regions;

[0123] The rate of change of albedo of each sampled region in the current image frame is calculated by difference with the reference rate in the time window. In operation, first select the image frame to be analyzed, for example T4 frame, extract the pixel set of a certain region Z1 in it, then obtain the rate of change of the corresponding pixels of the region during T3-T5 frames, for example the rate of pixel (150, 150) in the red band is 0.0003 / day, and the average rate in the time window T2-T4 is 0.001 / day, then the difference is 0.0007 / day. The absolute value of the difference is judged, if the difference is greater than the set albedo jump threshold, the pixel point is removed. The jump threshold is set by referring to the standard deviation of the rate difference of all pixels in the entire image, for example the standard deviation of the rate difference of all pixels in a band is 0.0005, and the jump threshold can be set to 0.0006, which means that when the rate difference of a certain pixel point is greater than 0.0006, it is determined to be abnormal change. If the number of removed pixels in a region exceeds 30% of the total number of pixels in the region, the entire region is not retained, for example, Z2 region contains 100 pixel points, of which 40 pixel points have rate difference exceeding the threshold, then the entire Z2 region is excluded. Conversely, if only 12 pixels in region Z1 have rate difference exceeding the threshold, and the remaining 88 pixels remain stable, then Z1 is determined to be an albedo stable region. Finally, all sampled frames and their corresponding regions are traversed, and the region set that meets the above conditions is output, which is the albedo stable region.

[0124] S303: Call the albedo stable region, extract the spatial distribution information of the retained region, arrange the region segment data in time sequence, and generate a set of sampled retained region segments;

[0125] Read the two-dimensional coordinate information of each pixel point in the region one by one, aggregate the pixel points of the region determined to be stable in each frame into image segments, for example, in frame T3, region Z1 contains pixels (100, 101), (101, 101), (102, 101), and in frame T4, it contains (100, 101), (101, 101), (101, 102). Group and arrange them into two region segments according to time frame, while retaining time frame index and region number. On this basis, sort the segments in time sequence, that is, record the corresponding position of each region number in the image frame in the order of T1, T2, T3, …, Tn. If a region number remains stable in consecutive frames, it is recorded as a continuous stable segment, if it is interrupted, it is treated as a separate segment, for example, Z a and Z b correspond to frames T2-T4 and T6 respectively. When constructing the final segment set, all region segments are arranged into structured data items according to frame index, spatial position, pixel number, etc. For example, Z aCorresponding to 3 frames, containing 12 pixels, the center position is (100, 100), Z b Corresponding to 1 frame, containing 9 pixels, all area segment information is summarized in the above format to generate a sample reserved area segment set.

[0126] The specific steps of S4 are:

[0127] S401: Call the sample reserved area segment set, extract the solar incident angle and terrain slope value corresponding to the area, calculate the projection length in the solar direction based on the incident angle and slope value, identify the light projection range of the area, and judge whether there is a slope shielding condition to obtain the slope shielding judgment result;

[0128] The spatial coordinate range and corresponding time frame index of each segment area are extracted one by one. For each pixel position in the area, the solar incident angle and terrain slope value on the corresponding remote sensing image frame are obtained. The solar incident angle data can be calculated by the image shooting time and the solar azimuth angle, for example, the solar incident angle of the image taken at 10 o'clock on June 21, 2023 is 35 degrees. The terrain slope value is calculated by the terrain elevation model data, for example, the slope value of a certain area is 20 degrees. Then, based on the relative angle between the incident angle and the slope value, the projection length of the ground in the solar direction is calculated. Specifically, the shadow length on each pixel point is converted through the triangular relationship. Under the condition of 30-degree incident angle and 20-degree slope, the angle is 10 degrees. If the height difference of the area is 2 meters, the projection length is about 11.4 meters, which may correspond to 3 to 4 pixel lengths on the image. By matching the projection length with the pixel coverage of the original area, it is judged whether some pixels are shielded by the projection. If a pixel is in the projection area and its slope direction relative to the solar direction deviates from the incident angle direction, it is judged that slope shielding occurs. The shielding condition is set as follows: when the angle between the slope direction and the incident direction exceeds 90 degrees and falls within the projection coverage, it is recorded as a shielding pixel. If the pixel points (120, 120), (121, 120), and (122, 120) in a certain area segment meet the condition, it is considered that the area exists slope shielding. The judgment results of all area segments are summarized to output the slope shielding judgment results corresponding to each area.

[0129] S402: Based on the slope shielding judgment result, extract the reflection intensity variation amplitude of the area, and perform combined operation with the projection value to calculate the combined value of the area and obtain the shielding influence combined value;

[0130] Extract the reflection intensity value of each pixel in the multi-phase image frame in the shielding area, and respectively count the reflection intensity variation amplitude of these pixels. The variation amplitude calculation method is the maximum reflection intensity value minus the minimum value. For example, the reflection intensity of pixel (120, 120) in T1-T6 frames is 0.26, 0.28, 0.25, 0.29, 0.27, 0.26, and the variation amplitude is 0.29 minus 0.25, which is 0.04. Perform this statistical operation on all shielding pixels in each region to calculate the average reflection intensity variation amplitude of the region as a whole. For example, there are 25 pixels in region Z1, and the average variation amplitude is 0.045. Then combine this value with the projection length calculated in the previous step. The combination operation rule is the product of the two, which is used to measure the reflection change intensity amplification caused by shielding. For example, the average projection length of Z1 region is 4 pixels, and the variation amplitude is 0.045, so the combined value is 0.18. After completion, perform the same operation process on all regions to obtain their respective shielding impact combined values. Then bind the region number and the corresponding combined value to form a structured record for subsequent judgment.

[0131] S403: Call the shielding impact combined value, filter the regions whose combined value exceeds the shielding judgment threshold, extract the corresponding image position and organize the spatial information, and generate a shielding exclusion image position set;

[0132] Screen each region's combined value in turn. Set the shielding judgment threshold to divide the shielding impact significant and non-significant regions. This threshold is set according to the distribution of combined values in the whole image. For example, the maximum value of all region combined values is 0.24, the average value is 0.12, and the standard deviation is 0.05. Then set the threshold to 0.18 as the screening limit. Compare the combined value of each region with 0.18. If the combined value is higher than the threshold, the region is marked as a shielding impact significant region. Extract the spatial position information of the region on the image, including the starting coordinates, pixel range, and time frame index, and convert it into a coordinate structure that the image processing module can recognize. For example, record the shielding area of region Z2 in image T4 frame (110, 110) to (114, 114) as T4_110_110_114_114. Organize all image positions that meet the shielding conditions in turn to generate a shielding exclusion image position set. Each item in this set contains information such as image number, spatial range, and region number, which can be used by the image processing system for subsequent region shielding exclusion operations.

[0133] The specific steps of S5 are:

[0134] S501: Call the shadow exclusion image position set, mark the distribution state of the remaining area in the spatial grid, track the corresponding land boundary extension form in the image frame, analyze the spatial connection of the area according to the continuity characteristics of the boundary adjacent area, and judge the integrity and continuity of the boundary to obtain the boundary coherent area;

[0135] The remaining area in the image frame is subjected to spatial grid division processing, and the image is divided into a plurality of fixed size sub-blocks using an equal interval grid structure, for example, each sub-block is 40x40 pixels. Whether each sub-block contains a remaining area is marked. If there is a pixel point coordinate not in the shadow exclusion set, the corresponding sub-block is marked as "valid", otherwise it is marked as "invalid". After recording the spatial distribution state, the valid sub-blocks in each frame of image are subjected to land boundary identification, and the extension form of the land boundary is extracted. The boundary connection relationship between sub-blocks is traversed in an 8-neighborhood manner during extraction, a land contour graph is constructed, and whether the contour has abnormal conditions such as breakage, interruption, and abnormal extension is identified. For example, if the middle sub-block is invalid, it is considered that the boundary is broken. Then, based on the spatial continuity characteristics of the boundary adjacent area, the continuity is judged. The continuity judgment is based on the effective state of adjacent sub-blocks and the coherence degree of pixel structure. The continuity judgment standard is set as follows: if the number of valid states in adjacent 3 sub-blocks is ≥2, and the boundary pixel points have overlapping row or column pixel numbers of more than 50%, it is judged as spatial continuity. If the condition is continuously met in a certain area along the horizontal direction or vertical direction for more than 4 sub-blocks, it is considered that the area boundary has continuity. In this way, all the boundary structures formed by the remaining areas are analyzed frame by frame, and finally a set of areas meeting the boundary integrity and continuity conditions is formed. The output is the boundary coherent area.

[0136] S502: Based on the boundary coherent area, the area meeting the sampling structure continuity judgment condition is screened, the area meeting the condition is classified as a candidate block, and the sampleable candidate area is obtained.

[0137] The continuity of the internal pixel sampling structure of each region is checked. The continuity of the sampling structure is determined by analyzing the minimum distance between sampling points, the arrangement direction, and the coverage density. The continuity determination condition is set as follows: the distribution distance of the sampling points in the region along the main direction is less than or equal to 80 pixels, and the number of pixels covered by all the sampling points in the region is greater than 30% of the total number of pixels in the region. For example, in a region with a size of 3200 pixels (40x80 pixels), if the sampling points cover a total of 1100 pixels, the coverage rate is 34.4%, which meets the density condition. If these points are uniformly distributed along the horizontal direction and the adjacent distance is between 70 and 80 pixels, the arrangement distance condition is met. The region is classified as a sampling structure continuous region. All boundary coherent regions that meet the above conditions are extracted and summarized as a candidate block set. The candidate blocks are numbered and additional spatial coordinate index, frame index, and pixel coverage range information are attached. For example, the candidate block C01 is located in the T3 frame, the starting coordinates are (120, 80), the range is 40x80 pixels, the sampling point coverage ratio is 35%, and the sampling distance is 75 pixels. A structured table record is formed for all candidate blocks. The final output is a sampleable candidate region.

[0138] S503: Call the sampleable candidate region, organize the spatial index information, and match the corresponding image position to establish the sampleable cultivated land image block distribution information.

[0139] Read the frame index, coordinate starting point, pixel range, and number information of each region. Organize these information into a unified spatial index structure. Each record contains four tuples of image frame number, candidate block number, and spatial position (starting row, starting column, row range, and column range). For example, the spatial index of candidate region C02 in T4 frame is T4_C02_100_120_40_40. Then, match these spatial indexes with the position mapping of each frame in the original image data. Correspond the index information to the specific pixel block position in the image matrix. Mark the actual distribution position of each sampleable region in the image. Group the distribution results by image frame. Each frame of image forms a distribution record list. For example, the T3 frame corresponds to candidate blocks C01, C03, and C05, and the corresponding positions are (120, 80), (160, 80), and (200, 120), respectively. Finally, combine and organize the spatial distribution data of all sampleable candidate blocks in all image frames to establish the sampleable cultivated land image block distribution information, which is used for subsequent multi-frame data fusion or partition sampling scheduling.

[0140] Please refer to Figure 2 A remote sensing monitoring system for cultivated land protection, comprising:

[0141] The time sequence gray scale change extraction module obtains remote sensing image sequences of the cultivated land region at different time nodes, extracts gray scale values of the same block of land pixels and sorts them, calculates the gray scale difference between adjacent time nodes, judges the difference value in combination with the gray scale average value in the image frame, marks the pixels that exceed the gray scale fluctuation reference value in continuous multiple changes, and generates a pixel region set;

[0142] The high-frequency fluctuation region screening module calls the pixel region set, sets a frequency threshold and marks a high-frequency region, checks the spatial position and coverage relationship of the marked region, adjusts the sampling level of the dense sampling region, records the adjusted pixel sampling order, and generates a time sequence image sampling priority reference value;

[0143] The albedo fluctuation screening module calls the time sequence image sampling priority reference value, extracts the multi-band albedo sequence of the adjusted sampling region, calculates the albedo fluctuation rate in the time window, screens out the regions whose rate difference exceeds the albedo jump threshold, and obtains a sampling reserved region fragment set;

[0144] The slope shading exclusion module calls the sampling reserved region fragment set, extracts the solar incident angle and terrain slope value of the corresponding region, calculates the solar direction projection length, judges the shading condition, screens the regions that exceed the shading determination threshold in combination with the projection value and the region reflection intensity fluctuation amplitude, and obtains a shading excluded image position set;

[0145] The cultivated land image block extraction module calls the shading excluded image position set, marks the spatial distribution of the remaining region, tracks the land block boundary in the image frame, judges the boundary integrity and continuity, screens the regions that meet the sampling structure continuity, outputs the spatial index and position set, and generates the sampleable cultivated land image block distribution information.

[0146] The above is only a specific embodiment of the present application, but the protection scope of the present application is not limited thereto, any skilled person in the art can easily think of changes or replacements within the technical range disclosed by the present application, which should be covered within the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. A remote sensing monitoring method for farmland protection, characterized in that, Includes the following steps: S1: Obtain remote sensing image sequences of cultivated land areas at different time nodes, extract the pixel gray values ​​of the same plot and arrange them by time, calculate the ratio of the difference between adjacent frames to the mean gray value within the frame, mark the positions that continuously exceed the fluctuation reference ratio, and generate a pixel region set. S2: Call the pixel region set, set the interference frequency threshold to filter frequent regions, check the position and coverage relationship frame by frame, perform priority downgrading on pixels in densely sampled regions, and generate time-series image sampling priority reference values. S3: Call the time-series image sampling priority reference value, extract the multi-band albedo sequence of the corresponding region, set the time window and calculate the inter-frame variation rate, compare the current frame variation rate with the fluctuation value in the window, and obtain the set of sampled and retained region segments. S4: Call the sampled and retained area fragment set, extract the solar incidence angle and slope parameters, calculate the illumination projection length and combine it with the reflection intensity change range, filter the area with the combined value exceeding the occlusion threshold, and obtain the occlusion exclusion image location set; S5: Call the masking and exclusion image location set, mark the remaining area, track the extension status of the plot boundary, determine the integrity and continuity of the boundary, and generate the distribution information of sampleable farmland image blocks; The specific steps of S1 are as follows: S101: Obtain remote sensing image sequences of cultivated land areas at different time nodes, extract and sort the pixel gray values ​​of the same area, calculate the gray value difference between adjacent time nodes, and calculate the ratio of the gray value difference of pixels by combining the gray value mean within the image frame, and obtain the pixel gray value ratio sequence. S102: Based on the pixel grayscale ratio sequence, filter the pixel positions where the grayscale ratio changes exceed a set reference ratio at multiple consecutive time nodes, and collect them to obtain a set of grayscale fluctuation marker pixels; S103: Call the set of grayscale fluctuation marker pixels, analyze the spatial distribution characteristics of adjacent pixels, merge high fluctuation pixels, and form a set of pixel regions; The specific steps of S2 are as follows: S201: Call the pixel region set, set the pixel frequency threshold, mark the regions that exceed the threshold, check the spatial position of the marked regions frame by frame, analyze the coverage, set the threshold and filter the regions with higher frequency, and obtain the distribution of high-frequency marked regions. S202: Based on the distribution of the high-frequency labeled areas, perform a sampling downsampling operation at a selected location in the densely sampled area, while retaining the original sampling data level of the unlabeled area, and obtain the adjusted pixel sampling level; S203: Call the adjusted pixel sampling level, record the sampling sorting structure of the pixels, organize them in chronological order, and establish a time-series image sampling priority reference value.

2. The remote sensing monitoring method for farmland protection according to claim 1, characterized in that, The pixel region set includes pixel location, grayscale fluctuation ratio, and change continuity marker; the temporal image sampling priority reference value includes sampling region, sampling frequency threshold, and pixel sampling order; the sampled retained region fragment set includes albedo sequence, fluctuation rate, and jump threshold; the occlusion exclusion image location set includes solar incidence angle, terrain slope, projection length, and occlusion determination threshold; and the sampleable farmland image block distribution information includes spatial index, location set, and sampling coherent region.

3. The remote sensing monitoring method for farmland protection according to claim 1, characterized in that, The specific formula for calculating the grayscale difference between adjacent time nodes is as follows: ； in, This represents the grayscale difference of pixel (x, y) between the i-th and i+1-th frames of the image. This represents the grayscale value of pixel (x, y) in the i-th frame. This represents the average grayscale value of all pixels in the i-th frame of the image. This represents the total number of pixels in the image. It represents the algebraic sum of the grayscale value changes of all pixels in the image between the i-th and i+1-th frames.

4. The remote sensing monitoring method for farmland protection according to claim 1, characterized in that, The specific steps for S3 are as follows: S301: Call the time-series image sampling priority reference value, extract the multi-band albedo sequence corresponding to the adjusted sampling area, set the time window and calculate the albedo variation rate between consecutive frames, organize the variation rate of the time nodes, and establish the albedo variation rate sequence. S302: Based on the albedo change rate sequence, compare the difference between the change rate of the region in the current sampling frame and the fluctuation rate within the time window, and filter out regions where the difference exceeds the albedo jump threshold to obtain the albedo stable region. S303: Call the albedo stable region, extract the spatial distribution information of the retained region, organize the region fragment data according to the time series, and generate a sampled retained region fragment set.

5. The remote sensing monitoring method for farmland protection according to claim 4, characterized in that, The specific formula for calculating the rate of change of albedo is as follows: ； in, This represents the rate of change of albedo at pixel (x, y) in the i-th frame of the image on band b. This represents the albedo value of pixel (x, y) in band b of the i-th frame image. This represents the average albedo value of all pixels in band b of the image frame. This represents the albedo value of pixel (x, y) in band b within the selected reference time frame. This indicates the time interval between the current image frame and the reference image frame.

6. The remote sensing monitoring method for farmland protection according to claim 1, characterized in that, The specific steps of S4 are as follows: S401: Call the sampled and retained area fragment set, extract the solar incidence angle and terrain slope value of the corresponding area, calculate the projection length in the solar direction based on the incidence angle and slope value, identify the illumination projection range of the area, and determine whether there is slope shading conditions, and obtain the slope shading judgment result. S402: Based on the slope shading judgment result, extract the variation range of the reflection intensity of the area, and perform a combination operation with the projection value to calculate the combination value of the area and obtain the shading influence combination value. S403: Call the combined value of the occlusion effect, filter the regions whose combined value exceeds the occlusion judgment threshold, extract the corresponding image positions and organize the spatial information to generate a set of occlusion exclusion image positions.

7. The remote sensing monitoring method for farmland protection according to claim 1, characterized in that, The specific steps of S5 are as follows: S501: Call the masking exclusion image location set, mark the distribution status of the remaining area in the spatial grid, track the extension shape of the corresponding plot boundary in the image frame, analyze the spatial connection of the area based on the continuity characteristics of adjacent boundary areas, and judge the integrity and coherence of the boundary to obtain the boundary coherent area. S502: Based on the boundary continuity region, filter the regions that meet the continuous determination condition of the sampling structure, classify the regions that meet the condition into candidate blocks, and obtain the sampleable candidate regions. S503: Call the sampleable candidate region, organize the spatial index information, match the corresponding image position, and establish the distribution information of sampleable farmland image blocks.

8. A remote sensing monitoring system for farmland protection, characterized in that, A remote sensing monitoring method for farmland protection according to any one of claims 1-7, the system comprising: The temporal grayscale change extraction module acquires remote sensing image sequences of cultivated land areas at different time nodes, extracts and sorts the grayscale values ​​of pixels in the same plot, calculates the grayscale difference between adjacent time nodes, combines the grayscale mean within the image frame to make a ratio judgment on the difference, marks pixels that exceed the grayscale fluctuation reference value in multiple consecutive changes, and generates a pixel region set. The high-frequency variation region screening module calls the pixel region set, sets the frequency threshold and marks the high-frequency region, checks the spatial position and coverage relationship of the marked region, adjusts the sampling level of the densely sampled region, records the adjusted pixel sampling order, and generates a time-series image sampling priority reference value. The albedo fluctuation screening module calls the time-series image sampling priority reference value, extracts the multi-band albedo sequence of the adjusted sampling area, calculates the albedo change rate within the time window, filters out areas where the rate difference exceeds the albedo jump threshold, and obtains a set of sample retention area segments. The slope occlusion removal module calls the set of sampled and retained area segments, extracts the solar incidence angle and terrain slope value of the corresponding area, calculates the solar direction projection length, judges the occlusion conditions, and combines the variation of regional reflection intensity and projection value to filter areas that exceed the occlusion judgment threshold, thus obtaining the occlusion removal image location set. The farmland image block extraction module calls the masked and excluded image location set, marks the spatial distribution of the remaining area, tracks the block boundaries in the image frame, judges the integrity and continuity of the boundaries, filters the areas that meet the sampling structure continuity, outputs the spatial index and location set, and generates the sampleable farmland image block distribution information.

Images