A farmland surveying and mapping method based on remote sensing data

Abstract

The application relates to the technical field of remote sensing surveying and mapping, and specifically discloses a farmland surveying and mapping method based on remote sensing data, which comprises the following steps: acquiring multi-source remote sensing data and constructing a unified space-time base data set; calculating the time-phase reliability of each time phase, and screening out a boundary stable window set; constructing a field block boundary candidate field, a field ridge candidate field and a ditch candidate field; performing field block-field ridge-ditch co-constraint collaborative vectorization to obtain an initial surveying and mapping result; performing topological closure repair and geometric attribute assignment to obtain a topologically consistent surveying and mapping result; identifying a low-confidence area and calling an alternative time phase for back substitution correction to output final surveying and mapping results. Through time-phase reliability screening, three-element co-constraint vectorization and closed-loop back substitution correction, the detection precision and topological consistency of farmland boundaries under complex terrains are significantly improved, the workload of manual review is reduced, standardized results that can be directly used in engineering are output, and the method can be widely applied to the fields of high-standard farmland construction acceptance, irrigation and drainage engineering configuration and the like.

Description

Technical Field

[0001] This invention relates to the field of remote sensing mapping technology, and more specifically to a method for surveying farmland based on remote sensing data. Background Technology

[0002] In agricultural modernization, high-standard farmland construction, and arable land protection, obtaining accurate and topologically consistent information on farmland elements such as field boundaries, ridges, and ditches is crucial. Traditional field measurement methods are highly accurate but inefficient and costly, making them difficult to implement on a large scale. With the development of remote sensing technology, using remote sensing imagery for farmland mapping has become the mainstream trend.

[0003] Currently, the mainstream technical approach for farmland mapping based on remote sensing data is as follows: First, acquire high-resolution optical remote sensing images from a single period or multiple time phases, or perform a simple fusion of optical images with synthetic aperture radar (SAR) images. Then, use image classification, semantic segmentation, or edge detection algorithms to identify farmland areas or field boundaries. Finally, through post-processing steps such as vectorization, edge smoothing, morphological restoration, and manual correction, generate the final field outline vector map. Some improved schemes introduce digital elevation models (DEMs) or phenological information to help improve identification accuracy.

[0004] However, existing technologies have significant limitations when applied to areas with complex terrain and fragmented plots, such as paddy fields in hilly or mountainous areas. On the one hand, these areas are often cloudy and foggy, making it difficult to acquire high-quality optical imagery. Furthermore, the reflectance characteristics of paddy fields change drastically during different phenological stages, such as irrigation, transplanting, and canopy closure, making it difficult to reliably extract the true geometric boundaries of the fields from single or simply fused remote sensing data. On the other hand, existing methods typically identify and extract fields, ridges, and ditches as independent objects, ignoring the inherent spatial topological constraints between them. This approach easily leads to problems such as unclosed field boundaries, adjacent fields adhering or incorrectly segmented, and broken or mismatched linear features of ridges and ditches with field boundaries. The resulting mapping results cannot meet the engineering application requirements for high-standard farmland acceptance and irrigation / drainage engineering design, which demand high data geometric accuracy and topological consistency. Summary of the Invention

[0005] The purpose of this invention is to provide a farmland mapping method based on remote sensing data, so as to solve the technical problems of unstable boundaries and inconsistent topological relationships in existing remote sensing farmland mapping methods under complex terrain conditions, especially in fragmented paddy field areas in hilly and mountainous areas.

[0006] The objective of this invention can be achieved through the following technical solution: a farmland mapping method based on remote sensing data, comprising the following steps:

[0007] Acquire multi-source remote sensing data and construct a unified spatiotemporal basis dataset;

[0008] The unified spatiotemporal basis dataset is calculated and filtered to obtain a set of boundary stability windows;

[0009] The boundary stable window set is subjected to multi-source fusion processing to obtain a fusion candidate field;

[0010] Perform co-constraint vectorization on the fusion candidate fields to obtain the initial farmland mapping results;

[0011] Based on the initial farmland mapping results, topology closure repair and geometric attribute assignment are performed to obtain topologically consistent mapping results;

[0012] The topologically consistent mapping results are identified to obtain low-confidence regions; and alternative time phase windows are called to perform back-substitution correction on the low-confidence regions to obtain corrected farmland mapping results.

[0013] Based on the corrected farmland survey results, the final farmland survey results are generated.

[0014] As a preferred embodiment of the present invention, acquiring multi-source remote sensing data and constructing a unified spatiotemporal basis dataset includes:

[0015] Acquire multi-temporal optical remote sensing data, multi-temporal synthetic aperture radar remote sensing data, high-resolution remote sensing data, and elevation remote sensing data of the target farmland area;

[0016] Preprocessing is performed on multi-temporal optical remote sensing data, multi-temporal synthetic aperture radar remote sensing data, high-resolution remote sensing data, and elevation remote sensing data respectively.

[0017] The preprocessed multi-source remote sensing data were standardized using a unified spatiotemporal reference.

[0018] Multi-source remote sensing data that have undergone unified spatiotemporal benchmark standardization processing are organized layer by layer to form a unified spatiotemporal base dataset.

[0019] As a preferred embodiment of the present invention, the unified spatiotemporal basis dataset is calculated and filtered to obtain a boundary stability window set, including:

[0020] For each time phase in the unified spatiotemporal basis dataset, calculate the time phase reliability separately;

[0021] Based on the phase reliability of each phase, the phase reliability of all phases are arranged sequentially in chronological order to generate a phase reliability sequence;

[0022] Based on the temporal reliability sequence, all temporal phases are sorted from high to low reliability, and the K temporal phases with the highest reliability are selected as the boundary stability window set, where K is the preset number of windows, and the value of K ranges from 3 to 5.

[0023] As a preferred embodiment of the present invention, the boundary stable window set is subjected to multi-source fusion processing to obtain a fusion candidate field, including:

[0024] For each time phase within the boundary stability window set, high-resolution optical remote sensing images are acquired and analyzed to obtain the optical image edge response, scattering discontinuity response, and elevation micro-topography edge response.

[0025] The edge response of optical images, the scattering discontinuity response, and the elevation micro-topography edge response are weighted and fused to generate candidate fields for field boundaries.

[0026] For each time phase within the boundary stability window set, high-resolution optical remote sensing images are acquired and analyzed to obtain the first optical linear enhancement response and elevation bulge response.

[0027] The first optical linear enhancement response and the elevation bulge response are weighted and fused to generate candidate field ridges;

[0028] For each time phase within the boundary stability window set, high-resolution optical remote sensing images are acquired and inverted to obtain the inverted images; the inverted images are analyzed to obtain the second optical linear enhancement response and elevation depression response.

[0029] The second optical linear enhancement response and the elevation depression response are weighted and fused to generate a ditch candidate field;

[0030] Based on the candidate fields for field boundaries, field ridges, and ditches, a fusion candidate field is generated.

[0031] As a preferred embodiment of the present invention, co-constraint vectorization is performed on the fused candidate fields to obtain initial farmland mapping results, including:

[0032] Non-maximum suppression processing is applied to the candidate fields for field boundaries, field ridges, and ditches to obtain refined candidate fields;

[0033] Candidate nodes are extracted based on the refined candidate field to obtain the candidate node set;

[0034] Construct a candidate edge set based on the candidate node set;

[0035] Construct an undirected graph based on candidate vertex set and candidate edge set. Where V is the candidate node set, E is the candidate edge set, and the label set is defined. , Let L represent the elements that the edge belongs to: field boundary, field ridge, ditch, or no element. An energy function is defined to evaluate the quality of the label assignment scheme L. The energy function expression is: ,in, For data items, These are constraints, including field boundary priority constraints, shared boundary uniqueness constraints, field ridge continuity constraints, ditch connectivity constraints, and field closure constraints. For the set of constraints;

[0036] The optimal label assignment result for the candidate edges is obtained by minimizing the energy function using the graph cut algorithm.

[0037] Based on the optimal label assignment results of the candidate edges, candidate edges assigned as field boundary labels are extracted to generate field polygons, candidate edges assigned as field ridge labels are extracted to generate field ridge lines, candidate edges assigned as ditch labels are extracted to generate ditch lines, and a topological relationship table between the elements is constructed to output the initial farmland survey results.

[0038] As a preferred embodiment of the present invention, candidate nodes are extracted based on a refined candidate field to obtain a candidate node set, including:

[0039] A confidence threshold is set, and pixel locations with a confidence level greater than the threshold are extracted from the refined candidate fields for field boundaries, field ridges, and ditches, serving as the first candidate node set. An 8-neighborhood connectivity analysis is performed on the initial candidate node set, and a breadth-first search algorithm is used to cluster interconnected nodes into candidate line segments. The endpoints of each candidate line segment are extracted and added to the second candidate node set. Intersections of three or more skeleton pixels within an 8-neighborhood on the skeleton line are identified and added to the third candidate node set. The first, second, and third candidate node sets are then merged to obtain the final candidate node set.

[0040] As a preferred embodiment of the present invention, the step of constructing a candidate edge set based on a candidate node set includes:

[0041] For every two candidate nodes in the candidate node set, calculate the spatial distance. If the spatial distance is less than the preset connection radius, generate a connection path using the Bresenham straight line algorithm and calculate the average confidence of the pixels on the path in the field boundary candidate field, field ridge candidate field, and ditch candidate field. If the maximum value of the average confidence of all pixels on the path is greater than the preset edge confidence threshold, establish a candidate edge between the two candidate nodes and record the starting point, ending point, path pixel coordinate sequence, and average confidence of each candidate field to generate a candidate edge set.

[0042] As a preferred embodiment of the present invention, the step of minimizing the energy function using a graph cut algorithm to obtain the optimal label assignment result for the candidate edges includes:

[0043] Construct a directed graph containing the source node, sink node, and all candidate edge nodes. Add directed edges from the source node and to the sink node to each candidate edge node. The capacity of each edge node corresponds to the cost of assigning field boundary labels and empty labels, respectively. Add bidirectional edges to candidate edge node pairs with constrained associations. The capacity of each bidirectional edge node corresponds to the penalty cost of inconsistent labels.

[0044] The minimum cut of the directed graph is solved using the maximum flow-minimum cut algorithm, dividing the graph into source and sink sides. Based on the minimum cut result, candidate edges located on the source side are assigned field boundary labels, and candidate edges located on the sink side are assigned empty labels. For candidate edges already assigned as field boundary labels, a secondary determination is made based on their relative confidence in the field ridge candidate field and the ditch candidate field. If the confidence in the field ridge is greater than that in the ditch candidate field, it is assigned as a field ridge label; otherwise, it is assigned as a ditch label. The optimal label assignment result for the candidate edges is obtained accordingly.

[0045] As a preferred embodiment of the present invention, topological closure repair and geometric attribute assignment are performed based on the initial farmland mapping results to obtain topologically consistent mapping results, including:

[0046] Topological defects are detected and located in the initial farmland mapping results, and a list of topological defects is generated.

[0047] Based on the list of topological defects, topological closure repair and regularization are performed on the field polygons, field ridge lines and ditch lines to obtain the repaired field polygons, field ridge lines and ditch lines;

[0048] Geometric attributes are calculated and assigned to the repaired field polygons, field ridge lines, and ditch lines to obtain the repaired and assigned field vectors, field ridge line vectors, and ditch line vectors.

[0049] Construct a topology table of field-field ridge-ditch relationship, integrate the repaired and assigned field vectors, field ridge line vectors, ditch line vectors and topology table, and output topology-consistent mapping results.

[0050] As a preferred embodiment of the present invention, the step of identifying the topologically consistent mapping results to obtain low-confidence regions, and then calling alternative time-phase windows to perform back-substitution correction on the low-confidence regions to obtain corrected farmland mapping results includes:

[0051] Based on the preset confidence assessment rules, low-confidence regions in the topological consistency mapping results are identified. The low-confidence regions include boundary offset regions, area anomaly regions, and linear facility breakage regions. The low-confidence regions are merged into low-confidence region patches.

[0052] For each low-confidence region patch, select 1-2 candidate time phases from the alternative time phases outside the boundary stability window set according to the time phase reliability and combine them with regional feature matching to form a back-substitution correction time phase set.

[0053] Using remote sensing data from the back-substitution correction time-phase set, candidate fields for field boundaries, field ridges, and ditches are reconstructed within low-confidence patch areas. Existing vector features at patch boundaries are used as fixed constraints to re-perform co-constraint vectorization and generate local correction results.

[0054] The local correction results are merged with the topology-consistent mapping results to smooth the transition at the boundary junctions, check and repair topology conflicts, and update the confidence level.

[0055] Repeat the above steps for iterative optimization until the total area of ​​newly identified low-confidence areas is less than the preset area threshold or the preset maximum number of iterations is reached, and generate the corrected farmland mapping results.

[0056] The beneficial effects of this invention are:

[0057] This invention constructs a temporal reliability model, integrating cloud cover ratio, phenological separability score, multi-source image geometric consistency score, and linear facility saliency score to select a set of boundary stability windows that are truly suitable for geometric boundary measurement. This ensures the stability of boundary extraction from the data source and effectively avoids boundary oscillation, closure failure, and area deviation problems caused by clouds, phenological changes, poor lighting conditions, etc. It is especially suitable for fragmented paddy field areas such as hilly and mountainous areas.

[0058] This invention constructs candidate fields for field boundaries by fusing optical image edge responses, scattering discontinuity responses, and elevation micro-topography edge responses; and constructs candidate fields for field ridges and ditches by combining optical linear enhancement responses with elevation convex and concave responses. This multi-source fusion strategy can integrate complementary information from different time phases, different sensors, and different land cover features, significantly improving the ability to extract boundaries and linear features under complex conditions such as crop shading, shadows, and changes in illumination.

[0059] This invention constructs a graph optimization model and applies field boundary priority constraints, shared boundary uniqueness constraints, field ridge continuity constraints, ditch connectivity constraints, and field closure constraints. By minimizing the energy function, it obtains the globally optimal label assignment for each edge and finally outputs closed field polygons, continuous field ridges and ditch lines, and generates a field-field ridge-ditch topology table. This ensures the topological consistency and engineering usability of the initial mapping results, lays a high-quality data foundation for subsequent topology repair and back-substitution correction, and significantly improves the automation level and reliability of farmland mapping in complex terrain.

[0060] This invention detects and repairs topological defects in initial mapping results, including forcibly sealing unclosed boundaries, unifying shared boundaries, merging fragmented areas, and correcting broken connections between field ridges and ditches, ensuring that the output meets engineering topology requirements. Simultaneously, it calculates the average area, perimeter, centroid, slope, and elevation of field plots, as well as the average width, standard deviation, and coefficient of variation of field ridges and ditches, enhancing the completeness and practical value of the results.

[0061] This invention identifies low-confidence areas through rules such as boundary offset, area anomalies, and linear facility breaks. It adaptively selects the optimal time phase from candidate time phases for local re-extraction and correction. Combined with fixed boundary constraints and a graph optimization model, a closed-loop optimization process is formed, significantly improving adaptability to complex local areas, reducing manual editing costs, and enhancing the reliability and consistency of the overall surveying results. The final surveying results include topologically consistent field polygons, field ridge lines, and ditch lines vector layers, as well as a complete field-ridge-ditch topological relationship table, assigned standardized attribute fields and confidence levels. The results can be directly imported into geographic information systems for high-standard farmland construction, arable land protection, and irrigation and drainage planning, supporting practical business operations such as project acceptance, area calculation, and design applications. Attached Figure Description

[0062] The invention will now be further described with reference to the accompanying drawings.

[0063] Figure 1 This is a flowchart of the farmland surveying method of the present invention.

[0064] Figure 2 This is a comparison diagram of the effects of the present invention and the prior art. Detailed Implementation

[0065] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0066] Please see Figure 1 As shown, this invention is a method for farmland mapping based on remote sensing data, comprising the following steps:

[0067] S1. Acquire multi-source remote sensing data and construct a unified spatiotemporal basis dataset;

[0068] Furthermore, the specific execution steps of S1 are as follows:

[0069] S101. Acquire multi-temporal optical remote sensing data, multi-temporal synthetic aperture radar remote sensing data, high-resolution remote sensing data, and elevation remote sensing data of the target farmland area;

[0070] Furthermore, multi-temporal optical remote sensing images are acquired from satellite remote sensing data distribution platforms or UAV aerial survey systems. These multi-temporal optical remote sensing data cover the entire rice growing season in the target farmland area, with a temporal resolution of 5-10 days and a spatial resolution of 10 meters. This includes Sentinel-2 or Landsat-8 / 9 satellite imagery. When acquiring images, priority is given to time phases with cloud cover below 20%, ensuring that images are included for at least the key phenological stages such as irrigation, transplanting, tillering, heading, and post-harvest bare soil stages.

[0071] Multi-temporal synthetic aperture radar (SAR) remote sensing images were acquired from the SAR satellite data distribution platform. The SAR remote sensing data and the SAR optical remote sensing data have the same time range, adopt an interferometric wide-swath mode, use a VV+VH dual polarization combination, have a spatial resolution of 10-20 meters, and a temporal resolution of 6-12 days, including Sentinel-1A / B satellite imagery.

[0072] High-resolution remote sensing imagery is acquired from commercial satellite data providers or UAV aerial photography systems. The high-resolution remote sensing imagery has a spatial resolution better than 1 meter, preferably sub-meter level satellite imagery of 0.3-0.8 meters, or a spatial resolution better than 0.2 meters, preferably 0.03-0.10 meters, from UAV orthophotos. When multiple high-resolution images cover the same area, data with a solar altitude angle greater than 40°, a shadow area ratio of less than 5%, and a stitching error between flight strips of less than one pixel is preferred.

[0073] Elevation remote sensing data is acquired from a geographic information data platform. The elevation remote sensing data has the same spatial extent as the high-resolution remote sensing image, including a digital surface model (DSM) or a digital elevation model (DEM). The ground sampling interval is preferably the same as or slightly coarser than that of the high-resolution remote sensing image, for example, 0.5-2 meters.

[0074] S102. Preprocess the multi-temporal optical remote sensing data, multi-temporal synthetic aperture radar remote sensing data, high-resolution remote sensing data and elevation remote sensing data respectively.

[0075] Furthermore, radiometric calibration, atmospheric correction, orthorectification, and cloud masking are sequentially performed on the multi-temporal optical remote sensing data. Atmospheric correction uses the 6S or MODTRAN model to convert image pixel values ​​into surface reflectance; orthorectification is based on a digital elevation model to eliminate geometric distortions caused by terrain; cloud masking uses Fmask or similar algorithms to generate cloud, cloud shadow, and snow / ice mask layers, and the cloud obstruction ratio for each temporal phase is statistically analyzed.

[0076] The multi-temporal synthetic aperture radar remote sensing data are sequentially subjected to radiometric calibration, multi-view processing, adaptive filtering for speckle removal, terrain correction, and geometric registration. Among them, terrain correction is based on a digital elevation model, using enhanced Lee filtering or Gamma-MAP filtering algorithms to suppress speckle, and projecting the backscattering coefficients from slant range geometry to geocoded ground geometry.

[0077] The high-resolution remote sensing image is subjected to radiometric and orthorectified correction to eliminate lens distortion and geometric deviations caused by terrain; the elevation remote sensing data is filtered and void-filled to eliminate the influence of non-ground features such as buildings and trees, and a digital elevation model representing the micro-topography of the surface is generated.

[0078] S103. Perform unified spatiotemporal benchmark standardization on the preprocessed multi-source remote sensing data;

[0079] Furthermore, using high-resolution remote sensing imagery as a benchmark, multi-temporal optical remote sensing data, multi-temporal synthetic aperture radar (SAR) remote sensing data, and elevation remote sensing data are uniformly converted into the same CGCS2000 Gauss-Kruger projection coordinate system as the high-resolution imagery. Using the boundary coordinates of the high-resolution imagery as a benchmark, vector clipping is performed on all multi-source remote sensing data to ensure complete spatial consistency across all data. Using the ground sampling interval (GSD) of the high-resolution imagery as a unified resolution, low-resolution optical, SAR, and elevation data are resampled, with bilinear interpolation used for reflectivity and scattering coefficient data, bilinear interpolation for elevation data, and nearest-neighbor interpolation for binary data such as cloud masks. Using the high-resolution imagery as a benchmark, sub-pixel-level geometric registration is performed on all multi-source remote sensing data, employing a polynomial transform model or triangulation correction to control the planar registration error to ≤0.5 pixels, ensuring a one-to-one correspondence between pixel center coordinates for all data. Multi-temporal optical cloud masks, invalid SAR regions, DEM anomalous elevation regions, and high-resolution imagery shadow regions are fused to generate a unified and effective mask.

[0080] S104. Organize the multi-source remote sensing data that have completed the unified spatiotemporal benchmark standardization process into layers to form a unified spatiotemporal base dataset.

[0081] Furthermore, the preprocessed and spatiotemporally unified multi-temporal optical time-series data, multi-temporal SAR time-series data, high-resolution remote sensing data, and elevation remote sensing data are organized into a unified spatiotemporal basis dataset D with a fixed structure. The dataset D includes a multi-temporal optical remote sensing data layer, a multi-temporal synthetic aperture radar remote sensing data layer, a high-resolution remote sensing image layer, an elevation remote sensing data layer, and a unified effective mask layer. All data layers meet the "six same" standard of being in the same coordinate system, the same projection zone, the same spatial range, the same pixel resolution, the same number of rows and columns, and the same pixel correspondence. It can be directly used for subsequent temporal reliability calculation and candidate field construction without additional geometric processing.

[0082] S2. Calculate and filter the unified spatiotemporal basis dataset to obtain a set of boundary stability windows;

[0083] Furthermore, the specific execution steps of step S2 are as follows:

[0084] S201. For each time phase t in the unified spatiotemporal basis dataset D, calculate the time phase reliability R. t Where t = 1, 2, ..., M, and M is the total number of time phase numbers;

[0085] Furthermore, the temporal reliability Rt is used to quantitatively evaluate the suitability of temporal data in geometric boundary measurement tasks. t The calculation formula is:

[0086] ;

[0087] Among them, C t To determine the cloud occlusion ratio of optical remote sensing data at different time phases, the proportion of pixels covered by clouds and cloud shadows to the total number of pixels is calculated from the cloud mask layer generated in the optical remote sensing data preprocessing stage of step S1.2. This yields C. t Value. C t The range of values ​​is C t The smaller the value, the richer the available optical information in the time phase; correspondingly, The term is used to characterize the proportion of effective observations in each time phase;

[0088] P t The phenological separability score is used to evaluate the spectral separability of field boundaries and background features at different time phases. The specific calculation method is as follows: First, the normalized water index and normalized vegetation index of the temporal optical remote sensing data are calculated separately; the formula for calculating the normalized water index is: The formula for calculating the normalized vegetation index is: ,in, These represent the reflectance in the green, red, and near-infrared bands, respectively. Then, based on a pre-defined field sample area, the difference between the normalized water index and the normalized vegetation index within the sample area is calculated, and this difference is normalized to... The interval is used to obtain the phenological separability score P. t Specifically, when the water index is high and the vegetation index is low during the irrigation period, the difference is relatively large, P t The value approaches 1; when the vegetation is fully covered during the canopy closure period, the difference is small, P t Approaching 0.

[0089] G t The geometric consistency score for multi-source images is used to evaluate the geometric registration consistency between multi-source images across different time periods. The specific calculation method is as follows: First, corner features are extracted from the temporal optical images, SAR images, and high-resolution images, for example, using the Harris corner detection algorithm or the SIFT feature point detection algorithm. Then, the extracted feature points from different images are matched, and the spatial distance error of the matched point pairs is calculated. Finally, the root mean square error (RMSE) of all matched point pairs is calculated, and the RMS error is normalized to a predetermined error threshold. For the interval, the normalization rule is: the smaller the root mean square error, the greater the G. t A higher value indicates better geometric consistency among multi-source images; conversely, a lower value indicates increased root mean square error (RMSE) if registration difficulties arise due to clouds, terrain, or registration algorithms. t The value decreased.

[0090] L t The linear feature saliency score is used to evaluate the identifiability of linear features such as field ridges and ditches in high-resolution remote sensing imagery over time. The specific calculation method is as follows: First, a linear filtering operator, such as Gabor filtering, is applied to the high-resolution remote sensing imagery to generate a linear feature response map. The Gabor filter kernel parameters are set based on the estimated field ridge width; for example, if the estimated field ridge width is Wr and the image ground sampling interval is GSD, then the filter kernel width is set to Wr / GSD pixels. Then, within a preset typical field sample area, the percentage of pixels in the linear feature response map whose response value exceeds a preset threshold is calculated. The preset threshold is, for example, 60% of the maximum response value. Finally, this percentage is used as the linear feature saliency score Lt. A higher percentage indicates clearer and more continuous linear features in the imagery. t The higher the value.

[0091] Let be the weighting coefficient, satisfying Furthermore, all weighting coefficients are non-negative real numbers. These weighting coefficients are obtained through sample calibration. The specific calibration method is as follows: 30-50 typical plots within the target farmland area are selected as calibration sample areas. The sample areas should cover different terrain types (flat land, gentle slopes, steep slopes), different plot shapes (regular plots, irregular plots), and different crop cover conditions. Professional technicians, based on prior knowledge, manually calibrate whether each time phase t is suitable for geometric boundary measurement, obtaining binary labels, where "1" indicates suitability and "0" indicates unsuitability.

[0092] The weight coefficients are determined using a grid search method. Specifically, candidate combinations of w1, w2, w3, and w4 are traversed within the interval [0, 1] with a step size of 0.1, ensuring that the sum of the four weights is 1. For each candidate weight, R is calculated for each calibration sample area and each time phase. t Value, and press R t Sort the values ​​from highest to lowest and filter out R. t The K time phases with the highest values ​​are used as the prediction window set. The proportion of time phases in the prediction window set that are manually labeled as "suitable for boundary measurement" is used as the basis for calculating the boundary F1 value. Simultaneously, based on the selected window set, subsequent steps S3-S5 are performed to calculate the field closure rate of the sample area output. The field closure rate refers to the proportion of successfully closed fields in the final generated fields. The weight combination that maximizes the weighted sum of the boundary F1 value and the field closure rate is selected as the final calibration result.

[0093] In one specific embodiment, the calibrated weighting coefficients are: w1=0.2, w2=0.4, w3=0.2, w4=0.2.

[0094] S202. Based on the time-phase reliability of each time phase, arrange the Rt values ​​of all time phases sequentially in chronological order to generate a time-phase reliability sequence. ;

[0095] It should be noted that the temporal reliability time series Each element corresponds to a comprehensive suitability score for an observation time phase, which intuitively reflects the quality of each time phase for farmland geometric mapping.

[0096] S203, Based on temporal reliability sequence According to the temporal reliability R t Sort all time phases from high to low, and select R. t The K phases with the highest values ​​are used as the boundary stability window set W, where K is the preset number of windows, and the value of K ranges from 3 to 5.

[0097] Furthermore, during the screening process, it is ensured that the boundary stability window set W contains at least one phase during the irrigation period or before and after transplanting, and one phase after harvest or during the bare soil period. Specifically, this is achieved by analyzing the phenological separability score P of each phase. t P was identified t The time phase when the value reaches a preset high threshold (e.g., greater than 0.7) is used as a representative of the irrigation and transplanting period to identify P. t Phases with values ​​below a preset low threshold (e.g., less than 0.3) are used as representatives of the harvest and bare soil periods. If the initial K selected phases do not simultaneously include both of the above two types of phases, then the candidate phases are selected according to R. t Values ​​are added sequentially in descending order until the phenological coverage constraint is met, thus avoiding phenological interference on the image extracted from the subsequent boundaries.

[0098] S3. Perform multi-source fusion processing on the boundary stable window set to obtain the fusion candidate field.

[0099] Furthermore, the specific execution steps of step S3 are as follows:

[0100] For each phase in the boundary stability window set W , To acquire its high-resolution optical remote sensing images The edge intensity map of the temporal image is extracted using the Canny edge detection operator. The high and low thresholds of the Canny operator are adaptively determined based on the image's grayscale distribution. Specifically, this is achieved by calculating the histogram of the gradient magnitude across the entire image and setting the low threshold... Set the high threshold to the 70th percentile of the histogram. Set as the 90th percentile of the histogram, the gradient magnitude is higher than Pixels are marked as strong edges, in and Pixels within the range are marked as weak edges, below Pixels in the weak edge are set to 0, pixels connected to the strong edge in the weak edge are retained, and the rest are discarded, resulting in the edge intensity map. .

[0101] Edge intensity maps were calculated for each of the K time phases in the boundary stability window set W. For each pixel x, take the maximum value among K time phases. The formula is: and will Linear normalization, the normalization formula is: The optical image edge response of the pixel is obtained. ,in, The full map Minimum and maximum values, optical image edge response Used to extract gradient features of field boundaries in high-resolution optical images.

[0102] It should be noted that this fusion strategy can integrate boundary information from different time periods and make up for the boundary loss caused by crop shading or poor light conditions in a single period. For example, the boundaries of paddy fields are clear during the irrigation period, while the boundaries are blurred due to crop cover during the canopy closure period. By taking the maximum value, the strong edges during the irrigation period can be preserved.

[0103] For each phase in the boundary stability window set W , SAR remote sensing images were acquired, and speckle suppression was performed using enhanced Lee filtering with a filter window size of 5×5 pixels and an edge preservation parameter of 0.5. The backscattering coefficients were estimated after filtering, and then the local gradient magnitude of the backscattering coefficients was calculated using the Sobel operator. The formula is: ,in, These are the horizontal and vertical gradients, obtained through convolution. The gradient magnitudes reflect abrupt changes in backscattering at points where water content or roughness changes.

[0104] Gradient magnitude maps were calculated for each of the K time phases in the boundary stability window set W. For each pixel x, take the maximum value among K time phases. The formula is: and will Linear normalization, the normalization formula is: The scattering discontinuity response of the pixel is obtained. ,in, The full map Minimum and maximum values, scattering discontinuous response Used to extract abrupt changes in backscattering in SAR images caused by variations in the water content or roughness of ground features. These abrupt changes are located at field boundaries.

[0105] Elevation remote sensing data is acquired from a boundary stable window set. The slope value of each pixel is calculated using the third-order inverse distance squared weighted difference method. The locations of abrupt slope changes are then extracted on the slope map using the Canny edge detection operator, resulting in a binary edge image. ,in, This indicates that the pixel is located at a point of abrupt change in slope, i.e., at the edge of a rapid change in terrain, such as a steep ridge, ridgeline, or gully edge. This indicates non-edge data. Since elevation data itself does not exhibit temporal variation, multi-temporal fusion is unnecessary. The edge response is directly normalized using binary linear normalization, with the normalization formula as follows: The elevation micro-topography edge response of the pixel is obtained. ,in, The full map Minimum and maximum values; elevation micro-topography response Used to extract edge features in elevation remote sensing data caused by micro-topographic undulations, which correspond to the location of field ridges (linear ridges) or ditches (linear depressions).

[0106] Optical image edge response Discontinuous scattering response Elevation micro-topography edge response Perform weighted fusion to generate candidate fields for field boundaries. The weighted fusion formula is: Where e represents the natural constant, For weight fusion.

[0107] In one specific embodiment, for areas with gentle terrain, optical information dominates, and the fusion weight calibration result is: This means it mainly relies on optical edge response; for mountainous areas, micro-topographic constraints are more important, and the fusion weight calibration result is... This means increasing the weight of elevation information to enhance micro-topographic constraints. It should be noted that the fusion weights can be pre-calibrated and stored according to different topographic zones, and the corresponding weight combination can be automatically selected according to the DEM slope statistics of the target area during actual processing.

[0108] For each temporal phase within the boundary stability window set, high-resolution optical remote sensing images are acquired, and linear structure enhancement is performed using the Frangi vascular enhancement filter. Specifically, the images are first Gaussian smoothed, then the Hessian matrix and eigenvalues ​​of each pixel are calculated, and the linear response value is calculated based on the eigenvalues. This filter highlights slender linear structures while suppressing patchy noise. To accommodate field ridges of varying widths, a multi-scale strategy is employed. A baseline scale is calculated based on the estimated ridge width and the image ground sampling interval. Three to five discrete scales are then uniformly selected within half to twice the baseline scale, and the linear response at each scale is calculated. The maximum value is taken as the response value for that pixel. After performing the above calculations for K temporal phases, for each pixel location, the maximum response value from all temporal phases is fused. The fused result is then linearly normalized to the 0-1 interval to obtain the first optical linear enhancement response.

[0109] Elevation remote sensing data was acquired and smoothed using median filtering. First, the direction of the field ridges was estimated. Within the calibrated sample area, the main direction of the ridges was statistically analyzed. If no clear main direction was found, a multi-directional scanning strategy was employed. For each pixel, the elevation difference response in four directions (0°, 45°, 90°, and 135°) was calculated, and the maximum value was taken as the response value. For each pixel, a certain distance was extended to both sides along a direction line perpendicular to the main direction of the ridge. This distance was equal to half the estimated ridge width divided by the ground sampling interval and rounded down. The elevation values ​​for the left, right, and center positions were obtained from the elevation remote sensing data using bilinear interpolation. The elevation difference of the bulge was calculated, which is the larger of the elevations on the left and right sides minus the center elevation. The minimum bulge threshold was set to 0.05 meters, and the maximum bulge threshold was set to 0.10 meters. If the elevation difference is less than the minimum threshold, the elevation elevation response of that pixel is 0; if the elevation difference is between the minimum and maximum thresholds, it is linearly mapped to the 0-1 range; if the elevation difference is greater than the maximum threshold, the response value is 1. If multi-directional scanning is used, the response value in all directions is calculated for each pixel, and the maximum value is taken as the final elevation elevation response.

[0110] The optical linear enhancement response and the elevation bulge response are weighted and fused to generate candidate field ridges. The weighted fusion formula is: ,in, To integrate weights, The first optical linear enhancement response and elevation bump response for each pixel x.

[0111] S3.3, For each phase in the boundary stability window set W To acquire its high-resolution optical remote sensing images Since ditches typically appear as dark linear features in optical images, such as water bodies and moist soil absorbing light, and Frangi filtering is sensitive to bright linear structures, the image is first inverted to convert the dark linear features into bright linear features. Then, multi-scale linear structure enhancement is performed on the inverted image using Frangi vascular enhancement filtering to obtain the response value of each pixel at the current time phase. The maximum value of the enhancement responses at K time phases is taken pixel by pixel, fused, and normalized to obtain the second optical linear enhancement response.

[0112] Elevation remote sensing data was acquired and smoothed using median filtering. First, the ditch direction was estimated. Within the calibration area, the main direction of the ditch was statistically analyzed. If no clear main direction was found, a multi-directional scanning strategy was employed. For each pixel, the elevation difference response in four directions (0°, 45°, 90°, and 135°) was calculated, and the maximum value was taken as the response value. For each pixel, a certain distance was extended to both sides along a line perpendicular to the main ditch direction. This distance was equal to half the estimated ditch width divided by the ground sampling interval and rounded. The elevation values ​​for the left, right, and center positions were obtained from the elevation remote sensing data using bilinear interpolation. The depression elevation difference was calculated by subtracting the center elevation from the larger of the left and right elevations. The minimum depression threshold was set to 0.05 meters, and the maximum depression threshold was set to 0.20 meters. If the indentation height difference is less than the minimum threshold, the elevation indentation response of that pixel is 0; if the indentation height difference is between the minimum and maximum thresholds, it is linearly mapped to the 0-1 range; if the indentation height difference is greater than the maximum threshold, the response value is 1. If multi-directional scanning is used, the response value in all directions is calculated for each pixel, and the maximum value is taken as the final elevation convexity response.

[0113] The optical linear enhancement response and the elevation depression response are weighted and fused to generate candidate fields for ditches. The weighted fusion formula is: ,in, To integrate weights, The second optical linear enhancement response and elevation bump response for each pixel x.

[0114] It should be further noted that the candidate fields for field boundaries, field ridges, and ditches are all raster images with the same coordinate system, resolution, projected coordinate system, and pixel alignment as the unified spatiotemporal basis dataset. The value of each pixel x is... This represents the confidence level that the pixel belongs to the field boundary, with values ​​normalized to [0,1], where 0 indicates that the pixel is definitely not a field boundary, and 1 indicates that the pixel is highly likely to be a field boundary; similarly, the value of each pixel x... This represents the confidence level that the pixel belongs to a linear feature along a field ridge; 0 indicates that the pixel is definitely not a field ridge, and 1 indicates that the pixel is highly likely to be a field ridge; the value of x for each pixel... The confidence level indicates whether the pixel belongs to a ditch or linear feature. 0 means that the pixel is definitely not a ditch, and 1 means that the pixel is very likely to be a ditch.

[0115] S4. Perform co-constraint vectorization on the fusion candidate fields to obtain the initial farmland mapping results;

[0116] Furthermore, the specific execution steps of S4 are as follows:

[0117] S401. Non-maximum suppression processing is applied to the candidate fields for field boundaries, field ridges, and ditches to obtain refined candidate fields;

[0118] Furthermore, non-maximum suppression (NMS) is performed on the three candidate fields to refine the candidate boundary responses and eliminate redundant response pixels. For each candidate field, all pixels are traversed, and the neighboring pixels along the gradient direction of each pixel are calculated. If the confidence value of the current pixel is less than the confidence value of its forward or backward neighboring pixels along the gradient direction, the confidence value of that pixel is set to zero. After NMS, the refined candidate field boundaries are obtained. Candidate fields along the ridges Candidate sites for ditches .

[0119] S402. Extract candidate nodes based on the refined candidate field to obtain a candidate node set;

[0120] Furthermore, set a confidence threshold. In one specific embodiment, From the refined candidate field , and Extracting those with confidence scores greater than a confidence threshold. The pixel position is used as the first candidate node set. For the initial candidate node set An 8-neighborhood connectivity analysis was performed, and a breadth-first search algorithm was used to cluster interconnected nodes into candidate line segments. The endpoints of each candidate line segment were then extracted and added to a second candidate node set. ; Identify intersections of 3 or more skeleton pixels within an 8-neighborhood on the skeleton line and add them to the third candidate node set. Finally, the candidate node set is obtained. ,and .

[0121] S403. Construct a candidate edge set based on the candidate node set;

[0122] Furthermore, for every two candidate nodes, the spatial Euclidean distance between them is calculated. If the distance is less than the preset connection radius... Then, the connection path between the two nodes is further examined. The value of the connection path is 5-10 times the ground sampling interval, preferably 8 times the ground sampling interval.

[0123] The Bresenham line algorithm is used to generate a path connecting two nodes. The positions of all pixels along the path are obtained, and the average confidence score of the pixels along the path in the candidate field is calculated. , and The average confidence level of all pixels on the path is determined if the maximum average confidence level of all pixels on the path is greater than the preset edge confidence threshold. Then, a candidate edge is established between the two nodes. The value range is 0.2-0.4, preferably 0.3, and attributes are recorded for each candidate edge e, including the starting node, ending node, path pixel coordinate sequence, and average confidence level on the candidate field of the field boundary. Average confidence level on candidate fields along the ridges Average confidence level on the candidate field of ditches .

[0124] S404. Constructing an undirected graph based on candidate vertex sets and candidate edge sets. Where V is the candidate node set, E is the candidate edge set, and the label set is defined. , Let L represent the elements that the edge belongs to: field boundary, field ridge, ditch, or no element. An energy function is defined to evaluate the quality of the label assignment scheme L. The energy function expression is: ,in, For data items, For constraint terms, For the set of constraints;

[0125] Data Items Used to measure the assignment of edge e as a label The lower the cost, the more reasonable the allocation. For each candidate edge e, the data item calculation formula is as follows, depending on its allocation label:

[0126] like ,but ;

[0127] like ,but ;

[0128] like ,but ;

[0129] like ,but ,in, Let be the penalty constant. The value range is 0.3-0.7, preferably 0.5;

[0130] The constraints include field boundary priority constraints, shared boundary unique constraints, field ridge continuity constraints, ditch connectivity constraints, and field closure constraints, which are used to penalize label assignments that violate the inherent geometric relationships between elements.

[0131] It should be noted that the definitions and specific calculation methods for each constraint term are as follows:

[0132] Field boundary optimization constraints refer to the requirement that field boundaries should be aligned with the spatial positions of field ridges or ditches. For fields assigned as... The confidence level of edge e in relation to candidate fields of ridges or ditches should not be too low; its penalty term ,in, The penalty coefficient is set to 1.0. The alignment threshold is set to 0.3.

[0133] The shared boundary uniqueness constraint means that the common boundary between adjacent fields can only be extracted once, and duplicate extraction is prohibited; its penalty term is... ,in, The penalty coefficient is set to 2.0. As an indicator function, if edge e has a shared conflict, that is, the same edge is marked as a boundary by two adjacent fields simultaneously, then ,otherwise .

[0134] The continuity constraint of field ridges means that the field ridges should form a continuous linear structure, and unnatural breaks are not allowed; its penalty term ,in, The penalty coefficient is set to 1.0. Let e ​​be the two endpoints of edge e. The degree of breakage at the endpoint is determined by the number of other connections at the endpoint. If the number of edges is 0, then If the value is 1 and the endpoint is a line feature endpoint, then... If the value is 2 and the endpoint is an interior point, then Other cases .

[0135] Canal connectivity constraints refer to the requirement that ditches maintain connectivity and form a linear network with flow direction; its penalty term ,in, The penalty coefficient is set to 1.0. For the degree of suspension, if the other connections at the endpoints are assigned as If the number of edges is 0, then If it is 1 and is the end of the network, then ,otherwise .

[0136] A field closure constraint means that each field polygon must be closed and its area must be greater than the smallest cartographic unit; its penalty term ,in, Let the closure penalty coefficient be 2.0. The area penalty coefficient is set to 1.0; This is a closure indicator function. If the distance between the first and last nodes of the polygon formed by several candidate edges is less than the closure tolerance, the value is 1; otherwise, it is 0. The closure tolerance is set to 1.5 times the ground sampling interval. The area of ​​the polygonal field plot. This is the area of ​​the smallest mapping unit, with a default value of 50 square meters.

[0137] S405. Minimize the energy function using the graph cut algorithm to obtain the optimal label assignment result for the candidate edges;

[0138] Furthermore, constructing directed graphs It includes the source node s, the sink node t, and all candidate node pairs for each candidate edge node. Add two directed edges (s,e) and (e,t), whose capacities correspond to the allocation of e as... and The cost is as follows: For each pair of candidate edge nodes (ei, ej), if there is a constraint associated, then a directed edge (ei, ej) and (ej, ei) are added, the capacity of which corresponds to the penalty cost for label inconsistency.

[0139] The maximum flow of a graph network H is solved using a maximum flow-minimum cut algorithm, such as the Boykov-Kolmogorov algorithm, to obtain the minimum cut. This algorithm iterates through growth, augmentation, and adoption phases to find augmenting paths from the source to the sink and pushes flow until no augmenting paths exist. The minimum cut divides the graph network into source and sink sides, and the cut edge set corresponds to the optimal label assignment. The final label of each candidate edge e is determined based on the minimum cut result. If e belongs to the source side, a label is assigned. If e belongs to the sink side, then assign a label. For those that need to be distinguished Based on the initial allocation, according to the situation, and A second determination is made based on the relative size; if Then allocate Conversely, allocation .

[0140] S406. Based on the optimal label assignment results of the candidate edges, extract the candidate edges assigned as field boundary labels to generate field polygons, extract the candidate edges assigned as field ridge labels to generate field ridge lines, extract the candidate edges assigned as ditch labels to generate ditch lines, construct a topological relationship table between elements, and output the initial farmland survey results.

[0141] Furthermore, all candidate edges assigned as field boundary labels are extracted to form a boundary edge set. Connectivity analysis is performed on the edges in the boundary edge set, and interconnected edges are clustered into boundary loops. For each boundary loop, if the distance between its first and last nodes is less than the closure tolerance, the first and last nodes are connected to form a closed polygon. If the area of ​​the polygon is greater than the area of ​​the smallest cartographic unit, it is output as a field polygon, and a unique identifier Field-ID is assigned to each field polygon.

[0142] Extract all candidate edges assigned as field ridge boundary labels to form a field ridge edge set. Perform connectivity analysis on the edges in the field ridge edge set, merge interconnected edges into field ridge line elements, extract their center lines, and assign a unique identifier Ridge-ID to each field ridge line.

[0143] Extract all candidate edges assigned as ditch boundary labels to form a ditch edge set. Perform connectivity analysis on the edges in the ditch edge set, merge interconnected edges into ditch line features, extract the center line, and assign a unique identifier Ditch-ID to each ditch line.

[0144] Identify adjacent field polygons, record shared edges and their lengths to form field-field topological relationships; identify the adjacency relationship between field ridge lines and field polygons, record the field ridge ID and the ID of the adjacent field, to form field-field ridge topological relationships; identify the adjacency relationship between ditch lines and field polygons, record the ditch ID and the ID of the adjacent field, to form field-ditch topological relationships.

[0145] The generated field boundary vector layers, field ridge line vector layers, ditch line vector layers, and topology table are stored in vector data format. The field boundary vector layers contain Field-IDs and polygon geometry information; the field ridge line vector layers contain Ridge-IDs and line geometry information; the ditch line vector layers contain Ditch-IDs and line geometry information; and the topology table stores the adjacency relationships between elements in a structured table format. Integrating the field boundary vector layers, field ridge line vector layers, ditch line vector layers, and topology table yields the initial farmland survey results.

[0146] S5. Based on the initial farmland mapping results, perform topology closure repair and geometric attribute assignment to obtain topology-consistent mapping results;

[0147] Furthermore, the specific execution steps of S5 are as follows:

[0148] S501. Detect and locate topological defects in the initial farmland survey results, and generate a list of topological defects;

[0149] Furthermore, the topological defects are detected element by element and segment by segment by traversing the field polygons, field ridge lines, and ditch lines in the initial farmland survey results. The detection contents include: field polygons that are not closed, self-intersecting, overlapping, or fragmented; field ridge lines that are broken, hanging, short-branched, or discontinuous; ditch lines that are isolated, have no flow direction, hanging endpoints, or are not connected; and adjacent field plots that share boundaries that are repeated, missing, or misaligned. The elements with the above-mentioned topological defects are marked one by one with their corresponding spatial locations to form a list of topological defects.

[0150] S502. Based on the list of topological defects, perform topological closure repair and regularization on the field polygons, field ridge lines and ditch lines to obtain the repaired field polygons, field ridge lines and ditch lines.

[0151] Furthermore, for unclosed field boundaries, within a closure tolerance range of 1.5-2 times GSD, the high-confidence paths along the candidate field Bp of the field boundary are automatically sealed to force the formation of closed polygons; for self-intersecting and overlapping field boundaries, the optimal path is retained according to the boundary confidence level, and redundant intersecting segments are eliminated; for fragments with an area smaller than the smallest cartographic unit, they are merged into adjacent fields according to spatial neighborhood and boundary similarity; for repeated or misaligned segments of shared boundaries, unified edge reduction processing is performed to ensure that the shared boundaries of adjacent fields are unique, consistent in location, and free from topological conflicts.

[0152] For breakpoints in the field ridge line where the break distance is less than 3 times the local average field ridge width, they are automatically connected along the high-confidence region of the field ridge candidate field Br to maintain the overall continuous linear shape of the field ridge; for suspended endpoints and isolated segments in the ditch line, they are extended along the ditch flow direction and the high-confidence path of the ditch candidate field Bd to connect to the main ditch network to ensure the connectivity of irrigation and drainage in the ditch; invalid short branches and isolated linear fragments with a length less than the threshold are removed so that the field ridges and ditches maintain a continuous, connected, and non-redundant engineered shape.

[0153] S503. Perform geometric attribute calculation and assignment on the repaired field polygons, field ridge lines and ditch lines to obtain the repaired and assigned field vectors, field ridge line vectors and ditch line vectors;

[0154] Furthermore, when calculating the geometric attributes of the restored field polygons, the geometric attributes include area, perimeter, centroid coordinates, minimum bounding rectangle, slope, and mean elevation. When calculating the geometric attributes of the restored field ridges and ditch lines, a skeleton extraction algorithm based on distance transformation or an active contour model is first used to extract and optimize the centerline: the linear features are rasterized, the distance from each raster to the boundary is calculated, local maxima are taken as skeleton points and connected to form the centerline, or the original line is adjusted to the position of the image gradient maxima using an active contour model; then sampling points are set along the centerline with a fixed step size of 0.5 meters, and a cross-sectional line perpendicular to the centerline is constructed at each sampling point, with the cross-sectional length being 3 times the estimated width of the feature; the grayscale profile of the high-resolution image and the elevation profile of the DEM are extracted from the cross-sectional line, using a... The Canny operator detects grayscale and slope abrupt changes, selecting the abrupt change locations as candidate points for the left and right boundaries. When the distance between the grayscale abrupt change and the elevation abrupt change location is less than 0.2 times the ground sampling interval, the location is determined as a boundary point. The spatial distance between the left and right boundary points is calculated as the cross-sectional width. When the local slope is greater than 5°, a slope correction coefficient is introduced. The specific correction formula is that the actual width equals the projected width divided by the cosine of the slope angle. The width values ​​of all cross-sections along the entire line are compiled into a sequence and smoothed using a median filter with a window size of 5. For outliers that deviate from the median width of adjacent cross-sections by more than 50%, they are replaced with the median width of adjacent cross-sections. Finally, the mean of the width sequence is recorded as the representative width, along with the minimum, maximum, and standard deviation, as well as the cumulative length and width variation coefficients along the centerline.

[0155] S504. Construct a topology table of field-field ridge-ditch relationship, integrate the repaired and assigned field vectors, field ridge line vectors, ditch line vectors and topology table, and output topology-consistent mapping results.

[0156] Furthermore, for each field ridge line, its spatial location and adjacency relationship with the field polygons are analyzed. Using a spatial overlay analysis method, the field ridge lines are buffered to both sides by 0.5 × Wridge, where Wridge is the width of the field ridge, generating buffer polygons. These buffer polygons are then overlaid with the field polygons; if the intersection area is greater than zero, an association is established. A topology table records the field ID, field ridge ID, and adjacency type, such as left-side adjacency or right-side adjacency, determined based on the relationship between the field ridge normal direction and the field center.

[0157] For each ditch line, an adjacency relationship with the field polygon is established using a method similar to that used for field ridges. The topology table records the field ID, ditch ID, and adjacency type, such as merging or adjacent. If a ditch line coincides with the field boundary and the ditch flow direction points inward into the field, it is marked as an irrigation inlet; if the flow direction points outward from the field, it is marked as a drainage outlet.

[0158] For each field polygon, analyze its spatial adjacency with surrounding field polygons. Using polygon overlay analysis: if two field polygons share a common boundary (i.e., the distance between boundary segments is less than 0.5 × GSD and the length is greater than 0), record the current field ID, adjacent field IDs, and the length of the shared boundary in the topology table. Also record the type of the shared boundary: field ridge boundary or natural boundary.

[0159] The topology-consistent mapping results include a topology-corrected field boundary vector layer containing geometric attributes such as area, perimeter, and centroid coordinates; a centerline-optimized field ridge line vector layer containing geometric attributes such as mean, minimum, maximum, and standard deviation of width; a centerline-optimized ditch line vector layer containing geometric attributes such as mean, minimum, maximum, and standard deviation of width; and a complete topology relationship table recording the relationships between fields and ridges, fields and ditches, and fields and adjacent fields.

[0160] It should be noted that step S5 is a post-processing of the initial farmland survey results output in step S4, which eliminates minor topological errors, assigns geometric attributes to fields, ridges and ditches, constructs topological relationships between elements, and outputs topologically consistent and attribute-complete survey results.

[0161] S6. Identify the topologically consistent mapping results to obtain low-confidence regions; and call the alternative time phase window to perform back-substitution correction on the low-confidence regions to obtain the corrected farmland mapping results.

[0162] Furthermore, the specific execution steps of S6 are as follows:

[0163] S601. Based on the preset confidence assessment rules, identify low-confidence areas in the topological consistency mapping results. The low-confidence areas include boundary offset areas, area anomaly areas, and linear facility breakage areas. Merge the low-confidence areas into low-confidence area patches.

[0164] Furthermore, the specific method for identifying low-confidence regions in the topologically consistent mapping results is as follows: For each boundary line segment in the field boundary vector, compare it with the field boundary candidate field Bp constructed in step S3. At each pixel position of the boundary line segment, calculate the offset distance between the actual position of the boundary line and the peak confidence position in the candidate field. If the offset distance exceeds a preset boundary offset threshold, the region where the pixel is located is marked as a boundary offset region. The boundary offset threshold is determined based on the ground sampling interval of the high-resolution remote sensing image, and its value ranges from 1.0 to 1.5 times the ground sampling interval. For example, when the ground sampling interval is 0.1 meters, the boundary offset threshold is set to 0.15 meters.

[0165] For each polygonal field, calculate its area A. currentSimultaneously, based on the boundary stability window set W selected in step S2, the reference value Areference for the area of ​​the field is estimated using multi-temporal data within the window set. Specifically, the estimation method is as follows: extract the boundary of the field from the remote sensing images of each temporal phase within the window set W, calculate the area at each temporal phase, and take the median or mean as the area reference value, using the formula... Calculate the area deviation rate ,like If the area deviation exceeds a preset threshold, the field is marked as an area with abnormal area. The area deviation threshold ranges from 3% to 5%.

[0166] For each field ridge line and ditch line, calculate its breakage rate. The breakage rate is defined as the ratio of the number of broken points on the line element to the theoretical number of continuous points. The specific calculation method is as follows: sample along the line element at a fixed step size, for example, 1 time interval of ground sampling, and count the proportion of sampling points whose confidence in the candidate field ridge Br or candidate field ditch Bd is lower than a preset threshold. The preset threshold is, for example, 0.3. This proportion is the breakage rate.

[0167] If the fracture rate exceeds a preset fracture rate threshold, the area where the line element is located is marked as a fracture area of ​​the linear facility. The fracture rate threshold ranges from 15% to 25%, preferably 20%.

[0168] The boundary offset regions, area anomaly regions, and linear facility breakage regions are spatially merged to form a set of low-confidence regions. Adjacent or overlapping low-confidence regions are merged to form continuous low-confidence region patches.

[0169] S602. For each low-confidence region patch, sort the candidate time phases outside the boundary stability window set according to the time phase reliability and combine them with regional feature matching to select 1-2 candidate time phases to form a back-substitution correction time phase set.

[0170] Furthermore, the temporal reliability sequence obtained from step S2 In the middle, select the one that is not included in the boundary stability window set W and has temporal reliability The higher phase is selected as the candidate phase, and R is obtained for all candidate phases. t Value, press R t The values ​​are sorted from high to low. For each low-confidence region, the main reasons for the decrease in confidence are analyzed. If the boundary of the region is clear during the irrigation period but the stable window does not include the phase, the phase near the irrigation period is selected first. If the boundary of the region is blurred due to crop shading, the phase after harvest or during the bare soil period is selected first. 1-2 alternative phases are selected for each low-confidence region patch to form a back-substitution correction phase set.

[0171] S603. Using remote sensing data from the back-substitution correction time-phase set, reconstruct candidate fields for field boundaries, field ridges, and ditches within the low-confidence area patches, and re-perform co-constraint vectorization using existing vector features at the patch boundaries as fixed constraints to generate local correction results.

[0172] Furthermore, for pixels within the low-confidence area patch, by back-substituting and correcting the resolution optical remote sensing images, SAR remote sensing images, and elevation remote sensing data in the time-phase set, and following the same execution process as step S3, local field boundary candidate fields, local field ridge candidate fields, and local ditch candidate fields are constructed.

[0173] For candidate fields within the low-confidence region patch area, the graph optimization collaborative vectorization process in step S4 is re-executed. During the local vectorization process, existing vector features at the boundary of the low-confidence region patch are used as fixed constraints to ensure smooth integration between the local correction results and the surrounding areas. Specifically, nodes on the low-confidence region boundary are set as fixed nodes, and their labels are inherited from the labels at the corresponding positions in the topologically consistent mapping results, and cannot be changed during the local optimization process.

[0174] By solving the local graph optimization problem, the local correction results within the low-confidence area patches are obtained, including the corrected field boundaries, field ridge lines, and ditch lines.

[0175] S604. Integrate the local correction results with the topology-consistent mapping results, smooth the transition at the boundary junctions, check and repair topology conflicts, and update the confidence level.

[0176] Furthermore, for the geometric features within low-confidence region patches, the original features in the topological consistency mapping results are replaced with features from the local correction results. At the boundaries of low-confidence region patches, if there are minor misalignments between the geometric features in the topological consistency mapping results and the local correction results, a smooth transition is achieved through linear interpolation or a fusion method based on elastic deformation to ensure geometric continuity at the boundaries. After fusion, a topological consistency check is performed on the low-confidence region patches and their neighborhoods. If new topological conflicts are found, such as newly added dangling lines or unclosed polygons, the topological repair method from step S5 is invoked for local correction. The confidence information of the fused region is updated. The confidence level of the region after back-substitution correction is marked as "corrected," and its confidence level is reassessed as A, B, or C based on the boundary offset, area deviation, and breakage rate after local correction.

[0177] S605. Repeat the above steps for iterative optimization until the total area of ​​newly identified low-confidence areas is less than the preset area threshold or the preset maximum number of iterations is reached, and generate the corrected farmland mapping results.

[0178] Furthermore, after completing one round of correction, step S601 is executed again to re-identify low-confidence regions. If low-confidence regions still exist and the maximum number of iterations has not been reached, the process returns to step S602, a candidate phase is selected again, and the next round of correction is performed. If the total area of ​​the newly identified low-confidence regions is less than the preset area threshold, such as 1% of the total surveyed area, or after the preset maximum number of iterations, such as 3, the iteration is stopped, and the corrected farmland survey results are output.

[0179] S7. Based on the corrected farmland survey results, generate the final farmland survey results.

[0180] Furthermore, the attribute fields of the corrected field vectors, field ridge vectors, and ditch vectors are standardized and assigned unique identifiers, areas, perimeters, lengths, average widths, maximum widths, minimum widths, width standard deviations, connectivity directions, and confidence levels. The topology table is organized into a relational data table, recording the association identifiers between fields, between fields and field ridges, and between fields and ditches, as well as the source field ID, target object type, target object ID, relationship type, and contact length.

[0181] Based on boundary offset, area deviation rate, fracture rate and width variation coefficient, the results are divided into three confidence levels: A, B and C, and a confidence level stratification map is generated.

[0182] The standardized surveying and mapping results are converted into the surveying and mapping database format, and metadata files containing project information, data sources, processing procedures, quality indicators and coordinate references are generated.

[0183] Perform geometric integrity, attribute integrity, topological integrity, and file integrity checks, and output the final farmland survey results.

[0184] Please see Figure 2As shown, compared with existing technologies, this invention demonstrates superior performance in terms of boundary detection accuracy, field closure rate, topological consistency, adaptability to complex terrain, and multi-element collaborative extraction accuracy. Specifically, this invention employs a temporal reliability model to screen a stable boundary window set, ensuring the stability of boundary extraction from the data source and effectively avoiding boundary oscillations and closure failures caused by phenological factors, cloud cover, and fog. Furthermore, it integrates optical, SAR, and elevation multi-source information to generate candidate fields for field boundaries, field ridges, and ditches. A graph optimization model is then used to apply common constraints such as field closure, boundary uniqueness, field ridge continuity, and ditch connectivity, achieving collaborative vectorization and topological consistency repair of the three elements. This solves problems such as boundary non-closure, breaks, and misalignments caused by element separation extraction in traditional methods. Further, through topological closure repair, automatic assignment of geometric attributes, and back-substitution correction in low-confidence regions, a closed-loop optimization is formed, significantly improving the surveying accuracy, automation level, and engineering adaptability under complex terrain conditions. The final output is a standardized surveying result with consistent topology and complete attributes, directly applicable to high-standard farmland acceptance and irrigation / drainage engineering design.

[0185] The above description is merely an example and illustration of the structure of the present invention. Those skilled in the art can make various modifications or additions to the specific embodiments described, or use similar methods to replace them, as long as they do not deviate from the structure of the invention or exceed the scope defined in the claims, all of which should fall within the protection scope of the present invention.

Claims

1. A method for farmland mapping based on remote sensing data, characterized in that, Includes the following steps: Acquire multi-source remote sensing data and construct a unified spatiotemporal basis dataset; The unified spatiotemporal basis dataset is calculated and filtered to obtain a set of boundary stability windows; The boundary stable window set is subjected to multi-source fusion processing to obtain a fusion candidate field; The boundary stable window set is subjected to multi-source fusion processing to obtain a fusion candidate field, including: For each time phase within the boundary stability window set, high-resolution optical remote sensing images are acquired and analyzed to obtain the optical image edge response, scattering discontinuity response, and elevation micro-topography edge response. The edge response of optical images, the scattering discontinuity response, and the elevation micro-topography edge response are weighted and fused to generate candidate fields for field boundaries. For each time phase within the boundary stability window set, high-resolution optical remote sensing images are acquired and analyzed to obtain the first optical linear enhancement response and elevation bulge response. The first optical linear enhancement response and the elevation bulge response are weighted and fused to generate candidate field ridges; For each time phase within the boundary stability window set, high-resolution optical remote sensing images are acquired and inverted to obtain the inverted images; the inverted images are analyzed to obtain the second optical linear enhancement response and elevation depression response. The second optical linear enhancement response and the elevation depression response are weighted and fused to generate a ditch candidate field; Based on the candidate fields for field boundaries, field ridges, and ditches, a fusion candidate field is generated; Perform co-constraint vectorization on the fusion candidate fields to obtain the initial farmland mapping results; Based on the initial farmland mapping results, topology closure repair and geometric attribute assignment are performed to obtain topologically consistent mapping results; The topologically consistent mapping results are identified to obtain low-confidence regions; and alternative time phase windows are called to perform back-substitution correction on the low-confidence regions to obtain corrected farmland mapping results. Based on the corrected farmland survey results, the final farmland survey results are generated.

2. The farmland surveying method according to claim 1, characterized in that, Acquire multi-source remote sensing data and construct a unified spatiotemporal basis dataset, including: Acquire multi-temporal optical remote sensing data, multi-temporal synthetic aperture radar remote sensing data, high-resolution remote sensing data, and elevation remote sensing data of the target farmland area; Preprocessing is performed on multi-temporal optical remote sensing data, multi-temporal synthetic aperture radar remote sensing data, high-resolution remote sensing data, and elevation remote sensing data respectively. The preprocessed multi-source remote sensing data were standardized using a unified spatiotemporal reference. Multi-source remote sensing data that have undergone unified spatiotemporal benchmark standardization processing are organized layer by layer to form a unified spatiotemporal base dataset.

3. The farmland surveying method according to claim 1, characterized in that, The unified spatiotemporal basis dataset is computed and filtered to obtain a set of boundary stability windows, including: For each time phase in the unified spatiotemporal basis dataset, calculate the time phase reliability separately; Based on the phase reliability of each phase, the phase reliability of all phases are arranged sequentially in chronological order to generate a phase reliability sequence; Based on the temporal reliability sequence, all temporal phases are sorted from high to low reliability, and the K temporal phases with the highest reliability are selected as the boundary stability window set, where K is the preset number of windows, and the value of K ranges from 3 to 5.

4. The farmland surveying method according to claim 1, characterized in that, Perform co-constraint vectorization on the fused candidate fields to obtain initial farmland mapping results, including: Non-maximum suppression processing is applied to the candidate fields for field boundaries, field ridges, and ditches to obtain refined candidate fields; Candidate nodes are extracted based on the refined candidate field to obtain the candidate node set; Construct a candidate edge set based on the candidate node set; Construct an undirected graph based on candidate node set and candidate edge set. Where V is the candidate node set, E is the candidate edge set, and the label set is defined. , Let L represent the elements that the edge belongs to: field boundary, field ridge, ditch, or no element. An energy function is defined to evaluate the quality of the label assignment scheme L. The energy function expression is: ,in, For data items, These are constraints, including field boundary priority constraints, shared boundary uniqueness constraints, field ridge continuity constraints, ditch connectivity constraints, and field closure constraints. For the set of constraints; The optimal label assignment result for the candidate edges is obtained by minimizing the energy function using the graph cut algorithm. Based on the optimal label assignment results of the candidate edges, candidate edges assigned as field boundary labels are extracted to generate field polygons, candidate edges assigned as field ridge labels are extracted to generate field ridge lines, candidate edges assigned as ditch labels are extracted to generate ditch lines, and a topological relationship table between the elements is constructed to output the initial farmland survey results.

5. The farmland surveying method according to claim 4, characterized in that, Based on the refinement of the candidate field, candidate nodes are extracted to obtain a candidate node set, including: A confidence threshold is set, and pixel locations with a confidence level greater than the threshold are extracted from the refined candidate fields for field boundaries, field ridges, and ditches, serving as the first candidate node set. An 8-neighborhood connectivity analysis is performed on the initial candidate node set, and a breadth-first search algorithm is used to cluster interconnected nodes into candidate line segments. The endpoints of each candidate line segment are extracted and added to the second candidate node set. Intersections of three or more skeleton pixels within an 8-neighborhood on the skeleton line are identified and added to the third candidate node set. The first, second, and third candidate node sets are then merged to obtain the final candidate node set.

6. The farmland surveying method according to claim 4, characterized in that, Construct a candidate edge set based on the candidate node set, including: For every two candidate nodes in the candidate node set, calculate the spatial distance. If the spatial distance is less than the preset connection radius, generate a connection path using the Bresenham straight line algorithm and calculate the average confidence of the pixels on the path in the field boundary candidate field, field ridge candidate field, and ditch candidate field. If the maximum value of the average confidence of all pixels on the path is greater than the preset edge confidence threshold, establish a candidate edge between the two candidate nodes and record the starting point, ending point, path pixel coordinate sequence, and average confidence of each candidate field to generate a candidate edge set.

7. The farmland surveying method according to claim 4, characterized in that, The optimal label assignment result for candidate edges is obtained by minimizing the energy function using the graph cut algorithm, including: Construct a directed graph containing the source node, sink node, and all candidate edge nodes. Add directed edges from the source node and to the sink node to each candidate edge node. The capacity of each edge node corresponds to the cost of assigning field boundary labels and empty labels, respectively. Add bidirectional edges to candidate edge node pairs with constrained associations. The capacity of each bidirectional edge node corresponds to the penalty cost of inconsistent labels. The minimum cut of the directed graph is solved using the maximum flow-minimum cut algorithm, dividing the graph into source and sink sides. Based on the minimum cut result, candidate edges located on the source side are assigned field boundary labels, and candidate edges located on the sink side are assigned empty labels. For candidate edges already assigned as field boundary labels, a secondary determination is made based on their relative confidence in the field ridge candidate field and the ditch candidate field. If the confidence in the field ridge is greater than that in the ditch, it is assigned as a field ridge label; otherwise, it is assigned as a ditch label. Based on this, the optimal label assignment result of the candidate edges is obtained.

8. The farmland surveying method according to claim 1, characterized in that, Based on the initial farmland mapping results, topology closure repair and geometric attribute assignment are performed to obtain topologically consistent mapping results, including: Topological defects are detected and located in the initial farmland mapping results, and a list of topological defects is generated. Based on the list of topological defects, topological closure repair and regularization are performed on the field polygons, field ridge lines and ditch lines to obtain the repaired field polygons, field ridge lines and ditch lines; Geometric attributes are calculated and assigned to the repaired field polygons, field ridge lines, and ditch lines to obtain the repaired and assigned field vectors, field ridge line vectors, and ditch line vectors. Construct a topology table of field-field ridge-ditch relationship, integrate the repaired and assigned field vectors, field ridge line vectors, ditch line vectors and topology table, and output topology-consistent mapping results.

9. The farmland surveying method according to claim 1, characterized in that, The low-confidence regions are identified by analyzing the topologically consistent mapping results. The alternative time phase window is then invoked to perform back-substitution correction on the low-confidence region, resulting in corrected farmland mapping results, including: Based on the preset confidence assessment rules, low-confidence regions in the topological consistency mapping results are identified. The low-confidence regions include boundary offset regions, area anomaly regions, and linear facility breakage regions. The low-confidence regions are merged into low-confidence region patches. For each low-confidence region patch, select 1-2 candidate time phases from the alternative time phases outside the boundary stability window set according to the time phase reliability and combine them with regional feature matching to form a back-substitution correction time phase set. By back-substituting and correcting the remote sensing data in the time-phase set, candidate fields for field boundaries, field ridges, and ditches are reconstructed within the low-confidence area patches. The existing vector features at the patch boundaries are used as fixed constraints, and co-constraint vectorization is re-executed to generate local correction results. The local correction results are merged with the topology-consistent mapping results to smooth the transition at the boundary junctions, check and repair topology conflicts, and update the confidence level. Repeat the above steps for iterative optimization until the total area of ​​newly identified low-confidence areas is less than the preset area threshold or the preset maximum number of iterations is reached, and generate the corrected farmland mapping results.

Images