A remote-sensing-based method for estimating the slope of a ploughed field

Abstract

This invention discloses a remote sensing-based method for estimating farmland slope, relating to the field of farmland slope technology. The method includes the following steps: acquiring image data and coarse-resolution digital elevation model (DEM) data; processing the image data to obtain a digital terrain model and orthophotos; visually interpreting the orthophotos to extract farmland distribution and shape indices; extracting topographic factors based on the DEM and filtering key factors using the variance inflation factor method; weighting the key factors using the entropy weighting method to obtain a comprehensive membership degree, combining a threshold to filter target slope points and overlaying farmland to obtain a preliminary farmland slope; performing a weighted average of the coarse-resolution DEM data to obtain the topographic slope; calculating the main slope aspect correction coefficient using the K-means clustering method on the DEM; fitting the preliminary farmland slope, topographic slope, main slope aspect correction coefficient, and shape index using a power function model to obtain a farmland slope estimation model, thus achieving the estimation of farmland slope.

Description

Technical Field

[0001] This invention relates to the field of field slope technology, specifically a method for estimating the slope of cultivated land based on remote sensing. Background Technology

[0002] Farmland slope is a core parameter in agricultural production, used for irrigation system design, soil and water conservation measures, and farmland quality assessment. It directly affects crop planting layout and field management efficiency. With the rapid development of remote sensing technology, estimating farmland slope at the regional scale based on remote sensing data has become an important alternative to traditional field measurements.

[0003] Currently, in the field of farmland slope estimation, the traditional method is based on a single coarse-resolution digital elevation model. This method relies on publicly available coarse-resolution digital elevation model data and combines remote sensing image interpretation with spatial overlay operations to obtain the farmland slope results.

[0004] However, the coarse-resolution digital elevation model relied upon by this traditional method is limited by spatial scale and cannot fully preserve micro-topographic features such as furrows and small depressions that directly affect the slope of the field. As a result, when calculating the slope of the cultivated land, it can only reflect the tilt trend of the macro-topography of the region and cannot accurately depict the slope details of the actual cultivated surface. Ultimately, the estimated results deviate significantly from the actual slope of the cultivated land and cannot meet the application requirements of high-precision slope data for refined agricultural management. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a remote sensing-based method for estimating the slope of cultivated land, thereby resolving the problems existing in the background technology.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a method for estimating the slope of cultivated land based on remote sensing, comprising the following steps:

[0007] Step S1: Obtain image data and coarse-resolution digital elevation model data for the area to be estimated;

[0008] Step S2: Rasterize and edit the image data to obtain a digital ground model and orthophoto map;

[0009] Step S3: Visually interpret the orthophoto map to obtain the distribution of cultivated land plots within the sample plots of the area to be estimated and calculate the shape index value of the cultivated land plots. Based on the digital terrain model, use the multi-window extraction method and slope formula to calculate and obtain multi-window terrain factors. Use the variance expansion factor method to screen the multi-window terrain factors and obtain key terrain factors.

[0010] Step S4: Based on the key terrain factors, the entropy weight method is used to perform weighted calculation to obtain the comprehensive membership degree of the multi-window terrain factors. The comprehensive membership degree of the multi-window terrain factors is compared with a preset threshold to select target slope points. The target slope points are spatially superimposed with cultivated land plots to obtain the preliminary field slope of the cultivated land plots.

[0011] Step S5: Based on the coarse-resolution digital elevation model data, the area-weighted average method is used to calculate the topographic slope of the cultivated land plot; the K-means clustering method is used to calculate the main slope aspect correction coefficient of the cultivated land plot using the digital ground model; the preliminary field slope, topographic slope, main slope aspect correction coefficient, and shape index value are fitted by regression analysis using a power function model to obtain the field slope estimation model;

[0012] Step S6: Input the shape index value, topographic slope and main slope aspect correction coefficient of the area to be estimated into the field slope estimation model, and output the preliminary field slope estimation results.

[0013] Preferably, the step of rasterizing and editing the image data to obtain a digital terrain model and orthophoto map includes the following specific steps:

[0014] The dense point cloud of the UAV was converted into a digital ground model with a resolution of 0.05m using ordinary Kriging interpolation. The interpolation formula is as follows:

[0015]

[0016] in, The estimated elevation values ​​for the raster to be interpolated. The number of known points in the interpolation neighborhood. For the first Interpolation weights for known points For the first Elevation values ​​of known points;

[0017] By combining UAV orthophotos, the elevation points of the digital ground model corresponding to the grids with gray values ​​<50 are removed; the grids with elevation values ​​>0.5m higher than the surrounding terrain are also removed. The minor undulations caused by median filtering are smoothed. The digital ground model is visually checked to ensure that 5cm-scale furrows and 20cm-scale depressions are retained without omission. Finally, the digital ground model and the corresponding orthophoto map are generated.

[0018] Preferably, the step of visually interpreting the orthophoto map to obtain the distribution of cultivated land plots within the sample plot of the area to be estimated includes the following steps:

[0019] Visual interpretation was performed based on orthophoto maps, and the following visual interpretation standards were set: cultivated land in black soil areas showed parallel ridge and furrow textures, cultivated land in plain areas showed uniform crop row textures, and the gray value range of cultivated land in orthophoto maps was 80-150. The boundaries of cultivated land plots were manually drawn in ArcGIS 10.8 software to generate initial cultivated land vector surfaces. A total of 12 cultivated land plots were extracted within the sample plots to ensure that there were no omissions or redundant drawings.

[0020] The initial boundary is optimized by first performing histogram equalization on the orthophoto image to enhance the grayscale difference between cultivated and non-cultivated land. Edge detection is then performed using the Canny operator, with the following operator parameters set: Gaussian filter standard deviation. =1.0, low threshold T1=50, high threshold T2=150;

[0021] The detected farmland edge vectors are superimposed with the visually interpreted initial boundaries, and then fused using the minimum distance method to finally generate the farmland plot boundary locations.

[0022] Preferably, the calculation of the shape index value of the cultivated land plot includes the following steps:

[0023] The formula for calculating the shape index of each arable land parcel is:

[0024]

[0025] in, This is a farmland plot shape index, reflecting the regularity of the plots. When =1, it is a square. The larger the plot of land, the more irregular its shape. The perimeter of the land parcel boundary. This refers to the area of ​​the land parcel.

[0026] Preferably, the step of calculating multi-window terrain factors based on the digital terrain model using a multi-window extraction method and slope formula includes the following steps:

[0027] The slope of each grid cell is calculated using a 3×3 moving window. The calculation formula is as follows:

[0028]

[0029] in, This represents the slope value for a single grid cell, in degrees. The elevation change rate in the x-direction is... The elevation change rate in the y-direction;

[0030] The terrain relief is the difference between the maximum and minimum elevation values ​​of the grid cells within the window; the surface cutting depth is the difference between the average and minimum elevation values ​​of the grid cells within the window. The formula for calculating the elevation variation coefficient is: ,in, The standard deviation of the elevation of the grid within the window. This represents the average elevation value of the grid within the window; A value > 0.01 indicates uneven terrain. A value < 0.005 indicates flat terrain;

[0031] To avoid missing micro-topographical scales with a single window, three windows with different pixel scales were selected to extract topographical factors: the window scales were 31×31 pixels, 41×41 pixels, and 51×51 pixels, corresponding to actual field ranges of 1.55m×1.55m, 2.05m×2.05m, and 2.55m×2.55m, respectively. The values ​​were based on experimental results, and the elevation variance converged within the 30-60 pixel window range. The 31-51 pixel window can cover the micro-topographical range of furrows and depressions within cultivated land, avoiding the omission of small-scale furrow factors by a single window.

[0032] For each farmland plot, five topographic factors are extracted under three windows to generate 15 sets of factor data.

[0033] Preferably, the step of filtering multi-window terrain factors using the variance inflation factor method to obtain key terrain factors includes the following steps:

[0034] The formula for calculating the variance inflation factor is:

[0035]

[0036] in, Let be the variance inflation factor of the i-th factor. denoted as the regression determination coefficient when the i-th factor is used as the dependent variable and all other factors are used as independent variables;

[0037] For each factor, a regression model is constructed with all other factors as independent variables. The model is then fitted using the least squares method to obtain the regression results for each factor. Substitute into the formula to calculate Factors with VIF > 10 are first removed; the VIF values ​​of the remaining factors are recalculated, and the factors with VIF < 10 are retained as key terrain factors, thus obtaining the key terrain factors.

[0038] Preferably, the step of using the entropy weighting method to perform weighted calculations based on the key terrain factors to obtain the comprehensive membership degree of the multi-window terrain factors includes the following specific steps:

[0039] For each window, a product-sum membership model is used. The formula for calculating the membership degree of a single window is as follows:

[0040]

[0041] in, Let k be the membership degree of the k-th pixel in each window. Let i be the weight of the i-th factor. This is the normalized value of the i-th factor in the k-th pixel;

[0042] The formula for calculating the comprehensive membership degree of multi-window terrain factors, which integrates the membership degrees of the three windows, is as follows:

[0043]

[0044] in, The multi-window terrain factor comprehensive membership degree of the k-th pixel. Let be the weight of the j-th window.

[0045] Preferably, the step of using K-means clustering to calculate the main slope aspect correction coefficient of the cultivated land plots includes the following specific steps:

[0046] The K-means clustering algorithm is used to aggregate discrete aspect values ​​into principal aspects, and the angular aspects are converted into unit circle coordinates to avoid the continuity problem between 0° and 360°. The calculation formula for aspect coordinate conversion is as follows:

[0047]

[0048] in, The slope angle, Unit circle coordinates;

[0049] Based on the slope aspect distribution characteristics, eight initial centers are set to correspond to eight directions: North, Northeast, East, Southeast, South, Southwest, West, and Northwest. The Euclidean distance from each slope aspect point to the eight centers is calculated, and the point is assigned to the nearest center. The center coordinates are updated, and this process is repeated until the center stabilizes. The number of iterations is 50, and the convergence threshold is 0.01. The formula for calculating the Euclidean distance is:

[0050]

[0051] in, and Let d be the coordinates of the d-th cluster center. and For the first The coordinates of the slope aspect point The Euclidean distance from each slope point to the 8 centers. The slope aspect point number;

[0052] Select the cluster center with the largest number of samples and convert its coordinates back to angular slope: , This reflects the angle between the long side of the cultivated land plot and the main slope direction. The smaller the angle, the greater the correction required for slope estimation. (Main slope direction correction factor) The calculation formula is: ,in, The angle of the longer side of the plot. The direction of the main slope of the quadrat.

[0053] Preferably, the step of performing regression analysis and fitting on the preliminary field slope, topographic slope, main slope aspect correction coefficient, and shape index value using a power function model to obtain the field slope estimation model includes the following specific steps:

[0054] By performing regression analysis on the field slope, topographic slope, main slope aspect correction coefficient, and shape index value using a power function model, a field slope estimation model is obtained. The calculation formula for the field slope estimation model is as follows:

[0055]

[0056] in, This represents the preliminary surface slope of a certain plot of farmland. For terrain slope, The shape index, The correction factor for the main slope aspect. and These are the model coefficients, obtained through parameter fitting.

[0057] Preferably, step S7 is also included, specifically:

[0058] Photon point cloud data of the area to be estimated is acquired. This data is then processed using a raster continuity-constrained filtering algorithm to obtain surface elevation control points for the area. An error correction model is constructed based on these control points. The preliminary field slope estimation results are input into the error correction model, and the resulting farmland slope estimation is output, thus achieving the estimation of farmland slope.

[0059] The photon point cloud data is processed using a filtering algorithm constrained by raster continuity to obtain the surface elevation control points of the area to be estimated. The initial raster size is calculated using the following formula:

[0060]

[0061] in, For grid width, This represents the highest elevation of the photon point within the sample plot area. This represents the lowest elevation of the photon point within the sample plot area;

[0062] Using the 8-connectivity standard, the connected components of all potential signal grids are labeled. The formula for calculating the grid continuity score is as follows:

[0063]

[0064] in, Raster continuity score; a higher score indicates stronger continuity between the raster and the surrounding terrain. The number of rasters in the connected region to the left of the raster to be detected. This represents the number of grid cells within the connected region to the right of the grid cell to be detected.

[0065] Estimated slope of the plot at each point The calculated slope of the field at 50 verification points was obtained through the preceding steps and was then measured using RTK (Real-Time Kinematics) to obtain the true values. Calculate the estimation error for each verification point. ,Think Independent variable, with With as the dependent variable, the formula for calculating the error correction model is:

[0066]

[0067] in, To estimate the error, Estimate the slope of the plot at each point. , and The correction coefficients are obtained through least squares regression fitting. =0.35, =-0.02, =0.005;

[0068] The error correction model is applied to the estimated slope of all plots, and the calculation formula is as follows:

[0069]

[0070] in, The results of the farmland surface slope estimation. To estimate the error.

[0071] This invention provides a method for estimating the slope of cultivated land based on remote sensing, involving machine learning and deep learning technologies, which has the following beneficial effects:

[0072] (1) This technical solution integrates multi-source data and constructs a high-precision digital ground model of the sample plot. It generates dense point clouds through UAV motion recovery structure and multi-view stereo vision technology, and combines photon point cloud preliminary denoising. It can completely preserve micro-topographic features such as ridges and depressions in cultivated land, and solves the problem that traditional single data sources are difficult to capture fine terrain and have insufficient accuracy, providing reliable basic data support for subsequent slope calculation.

[0073] (2) Based on the constructed high-precision digital ground model, the variance expansion factor method is used to quantify the multicollinearity of the topographic factors of cultivated land plots to eliminate redundancy. Then, the entropy weight method is used to objectively calculate the weight of key topographic factors, avoiding the problem of topographic factor interference caused by traditional subjective weighting. This makes the selection of topographic factors and weight allocation more in line with the actual topographic characteristics of cultivated land, accurately distinguishes the influence of target slope points and interference features, and further improves the accuracy of field slope calculation in sample plots.

[0074] (3) K-means clustering is introduced to determine the main slope aspect of cultivated land and calculate the correction coefficient. An improved power function model is constructed, which breaks through the limitations of the traditional model that only considers the topographic slope and shape index and ignores the topographic anisotropy. It can better adapt to the estimation of cultivated land slope in complex terrain areas, and make the model's slope prediction for different slope aspects more consistent with the actual tilt of the cultivated surface.

[0075] (4) Based on the constructed power function model, the surface elevation control points are obtained by using the continuity constraint filtering of photon point cloud grids, and a quadratic polynomial error correction model is constructed. This effectively eliminates the accumulation of systematic errors in large-scale estimation, solves the problem that the model is prone to accuracy decline due to data deviation when applied to the whole area, and further ensures the reliability of large-scale farmland slope estimation. Attached Figure Description

[0076] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0077] Figure 1 This is a flowchart illustrating the steps of a remote sensing-based method for estimating the slope of cultivated land proposed in this invention.

[0078] Figure 2 This is a step hierarchy diagram of obtaining the field slope estimation model in a remote sensing-based method for estimating farmland slope proposed in this invention.

[0079] Figure 3 This is a hierarchical diagram of the steps involved in obtaining the farmland slope estimation result in a remote sensing-based farmland slope estimation method proposed in this invention. Detailed Implementation

[0080] 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.

[0081] Please see Figures 1-3 The present invention provides a technical solution: a method for estimating the slope of cultivated land based on remote sensing.

[0082] S1: Obtain image data and coarse-resolution digital elevation model data for the area to be estimated.

[0083] This step requires acquiring image data and coarse-resolution digital elevation model data for the area to be estimated. The image data includes UAV aerial photographs and optical satellite imagery data. The specific data types and core parameters are as follows:

[0084] An industrial-grade quadcopter drone was selected, equipped with a mechanical global shutter aerial camera (model XCam1200s). Specific parameter settings and rationale are as follows: Spatial resolution is set to 5cm, which preserves micro-topographic features such as furrows and depressions in farmland; the minimum overlap rate for SfM (Structure of Motion) is 85%, ensuring that when SfM technology generates dense point clouds, the matching rate of corresponding points in adjacent drone aerial images is >90%, avoiding point cloud voids; the flight altitude is determined using a flight altitude formula to ensure that the ground resolution (GSD) meets the 5cm requirement. The calculation formula is as follows: ,in, The altitude of the drone relative to the ground. For the focal length of the aerial camera lens, For ground resolution, To determine the camera pixel size, RTK (Real-Time Kinematics) technology is used to eliminate flight errors, ensuring the final ground resolution remains at 5cm. Both lateral and forward overlap rates are set to 70% (usually lower than the minimum overlap rate, as actual flight errors exist, and the overlap rate design needs to reserve redundancy). The number of aerial photographs to be taken is calculated as: ceil(sample width ÷ (swath width × (1 - lateral overlap rate))) × ceil(sample length ÷ (swath width × (1 - forward overlap rate))), where ceil() is the floor function. For example, when the sample plot is 72m×72m, the XCam1200s camera has a 40m swath width at an altitude of 80m. At this time, the calculated number of aerial photos to be taken is 36. To avoid missing drone aerial photos, it is assumed that the probability of drone images failing to match due to differences in viewing angle is about 15%. Theoretically, with 36 aerial photos, the number of failed aerial photos is calculated as 36×15%=5.4. At least 6 redundant aerial photos need to be added to ensure that the number of effective aerial photos is ≥36. In reality, 42 drone aerial photos need to be taken.

[0085] SRTM DEM V3 data was selected as the coarse-resolution digital elevation model data, with a spatial resolution of 30m and an absolute elevation error of <16m (SRTM DEM official accuracy index). It was downloaded through the NASA EarthData platform to ensure coverage extending 1km beyond the boundary of the area to be estimated.

[0086] The Gaofen-7 satellite 1A-level dual-line array imagery with a spatial resolution of 0.65m was selected. The sample plots were typical representative sub-regions selected within the area to be estimated (covering the core terrain and cropping characteristics of the area to be estimated). Imagery with a time difference of less than 15 days from the UAV sample plot data was selected. Based on the area to be estimated and the width of the Gaofen-7 imagery, if the lateral length and flight path length of the area to be estimated are ≤20km and ≤20km (i.e., the coverage area is ≤20km×20km), one image can completely cover it. If the lateral length or flight path length of the area to be estimated is >20km or >20km (e.g., 10km×30km), multiple images are required for joint coverage. In this case, adjacent images are selected, requiring an inter-image overlap rate >10% (the overlap ratio of two adjacent images in spatial range) to ensure that subsequent image mosaicking is seamless.

[0087] The image data and the coordinates of the coarse-resolution DEM were uniformly converted to the WGS84 coordinate system (geodetic coordinate system). The specific processing is as follows:

[0088] The exterior orientation elements of the imagery are converted to WGS84 coordinates using POS data (GPS positioning accuracy ±0.1m) recorded by the UAV flight control system; ICESat-2ATL03 data uses the WGS84 coordinate system by default and does not require additional conversion; the coarse resolution DEM and satellite imagery are converted to the WGS84 coordinate system using the "Coordinate Conversion Tool" of ENVI5.6 software, with a conversion error of <0.5m.

[0089] S2: Rasterize and edit the image data to obtain a digital ground model and orthophoto map.

[0090] The image data is rasterized and edited to obtain a digital terrain model and orthophoto map, as specifically implemented below:

[0091] First, the SfM (Structure of Motion) + MVS (Multi-view Stereo Vision) technology workflow is used to convert the image data acquired in step 1 into a high-density, low-noise 3D point cloud. The specific implementation steps are as follows:

[0092] First, radiometric correction is performed on the 42 image data to eliminate grayscale unevenness caused by lens distortion and exposure differences. The formula for correcting camera distortion is:

[0093]

[0094] in, To correct the distance from the foreground image point to the principal image point, the distance is calculated using pixel coordinates from the UAV aerial photograph. To correct the distance from the image point to the principal point, and The radial distortion coefficient has the following values: , .

[0095] The Harris operator and the DoG operator are combined to extract feature points from UAV aerial images. The Harris operator detects corner points (inflection points of furrows in cultivated land and corner points of field boundaries), and the DoG operator detects edge points (edges of crop rows). The number of feature points extracted from a single UAV aerial image is ≥5000, and the matching rate of corresponding points between adjacent UAV aerial images is >90%.

[0096] Feature point matching employs Euclidean distance matching combined with RANSAC (random sampling consensus) algorithm to eliminate gross errors. The matching relationship is constrained by the fundamental matrix to ensure matching accuracy. Finally, the reprojection error of the fundamental matrix is ​​less than 1 pixel, ensuring the accuracy of subsequent SfM reconstruction.

[0097] Based on global SfM technology, the exterior orientation elements of UAV aerial photographs are solved and sparse point clouds are generated through the following steps:

[0098] Two UAV aerial images with the highest overlap (85%) were selected as seed images. The initial rotation matrix R and translation vector T were solved using fundamental matrix decomposition. The remaining 40 UAV aerial images were then added sequentially. The camera intrinsic parameters (focal length, principal point) and exterior orientation elements were optimized using bundle adjustment, requiring the reprojection error after adjustment to be <0.5 pixels. The object space coordinates (WGS84 coordinate system) of corresponding feature points were calculated using triangulation to generate a sparse point cloud with a point cloud density of approximately 5 points / m². 2 Key topographic features such as farmland boundaries and furrows should be preserved.

[0099] Based on sparse point clouds, a dense point cloud is generated using multi-view stereo vision technology. The specific process is as follows: For each UAV aerial image, a disparity space is constructed based on the sparse point cloud. The PatchMatch algorithm is used to estimate the depth value (unit: m) of each pixel, requiring a depth estimation error of <0.05m. Consistency checks are performed on the depth maps of all UAV aerial images (pixels with a depth deviation >0.1m are removed), and the images are fused to generate a dense point cloud. The formula for dense point cloud computing is: ,in, Point cloud density, in points / m 2 , This represents the total number of dense point clouds within the sample plot area. The area of ​​the sample plot is denoted as .

[0100] The dense point cloud from the drone is converted into a raster digital terrain model (DTM), and non-terrain interference is eliminated through editing. The specific implementation steps are as follows:

[0101] The spatial resolution of the DTM is set to 5cm. Ordinary Kriging interpolation is used to convert the dense point cloud of the UAV into a raster DTM with a resolution of 0.05m. The interpolation formula is:

[0102]

[0103] in, The estimated elevation values ​​for the raster to be interpolated. The number of known points in the interpolation neighborhood. For the first Interpolation weights for known points For the first Elevation values ​​of known points (dense point cloud of drones).

[0104] By combining UAV orthophotos (DOM), DTM elevation points corresponding to grids with gray values ​​<50 (tree shadow areas) are removed; grids with elevation values ​​>0.5m higher than the surrounding terrain are also removed. Median filtering (3×3 window) is used to smooth out minor undulations (elevation difference <3cm) caused by farmland stubble. Visual inspection is conducted to ensure that the DTM retains furrows at the 5cm scale and depressions at the 20cm scale, ensuring no omissions. Finally, the sample plot digital terrain model (DTM) and the corresponding orthophoto map (DOM) are generated.

[0105] S3: Visually interpret the orthophoto map to obtain the distribution of cultivated land plots within the sample plots of the area to be estimated and calculate the shape index value of the cultivated land plots. Based on the digital terrain model, use the multi-window extraction method and slope formula to calculate and obtain multi-window terrain factors. Use the variance expansion factor method to screen the multi-window terrain factors and obtain key terrain factors.

[0106] Farmland plots are the basic units for estimating field slope. It is necessary to first accurately separate farmland from non-farmland (such as roads, ditches, and areas shaded by trees) from the quadrat plots, and then calculate the plot morphology parameters. The specific implementation is as follows:

[0107] The orthophoto map DOM, generated synchronously with the sample plot DTM in step 2, is used. The DOM can intuitively reflect the texture characteristics of cultivated land. Visual interpretation is performed based on the orthophoto map, and the visual interpretation standards are set as follows: cultivated land in black soil areas shows parallel furrow texture (furrow spacing 0.5-1m), and cultivated land in plain areas shows uniform crop row texture. The gray value range of cultivated land in the DOM is 80-150. The boundaries of cultivated land plots are manually drawn in ArcGIS 10.8 software to generate the initial cultivated land vector surface, based on the upper and lower limits of the specified sample plot area (500m²). 2 -1000m 2 ) and the upper and lower limits of the area of ​​a single plot (50-200m²) 2 The theoretically reasonable range for estimating the number of plots within a sample plot (minimum theoretical number = lower limit of sample plot area) Maximum area of ​​a single plot, maximum theoretical number = maximum area of ​​sample plots The minimum theoretical number of plots is calculated to be 3, and the maximum theoretical number is 20. The theoretical reasonable range of the number of plots in the sample plot is (3-20). This implementation takes the middle value of the theoretical reasonable range and extracts a total of 12 cultivated land plots in the sample plot. It should be noted that 12 is not the only immutable number. Other numbers can be substituted under the following conditions: within the theoretical reasonable range (3-20) and the cultivated land integrity rate is >98%.

[0108] To address the issue of insufficient boundary smoothness in visual interpretation, the initial boundary is optimized. First, histogram equalization is applied to the DOM to enhance the grayscale difference between cultivated and non-cultivated land. Then, edge detection is performed using the Canny operator, with the following operator parameters set: Gaussian filter standard deviation. =1.0, low threshold T1=50, high threshold T2=150. This parameter combination can completely preserve the 5cm wide furrow boundary while removing noisy edges.

[0109] The detected farmland edge vectors are superimposed with the visually interpreted initial boundaries, and the minimum distance method is used for fusion (the operator detection results are retained when the boundary deviation is <0.1m). The final generated farmland plot boundary position error is <0.05m, which meets the spatial accuracy requirements for subsequent terrain factor extraction.

[0110] The formula for calculating the shape index of each arable land parcel is:

[0111]

[0112] in, This is a farmland plot shape index, reflecting the regularity of the plots. When =1, it is a square. The larger the plot of land, the more irregular its shape. The perimeter of the land parcel boundary. This refers to the area of ​​the land parcel.

[0113] Topographic factors are the core correlation indicators of field slope. Based on the aforementioned digital terrain model, multi-window extraction and slope formulas are used to calculate multi-window topographic factors, as detailed below:

[0114] Based on the influencing factors of farmland slope (elevation difference, local slope, and undulation), the following five initial factors were selected: Elevation: Original DTM raster elevation value, reflecting the elevation benchmark within the plot; Slope: Reflecting the local slope of the plot, related to the field slope; Topographic relief: Reflecting the local elevation difference within the plot, related to the influence of micro-topography (such as depressions and protrusions) on slope; Surface incision depth: Reflecting the depth of erosion of the plot, related to the furrow cutting characteristics of the black soil region; Elevation variation coefficient: Reflecting the dispersion of elevation data, related to topographic uniformity.

[0115] The slope of each grid cell is calculated using a 3×3 moving window. The calculation formula is as follows:

[0116]

[0117] in, This represents the slope value for a single grid cell, in degrees. The elevation change rate in the x-direction is... Let be the rate of change of elevation in the y-direction.

[0118] It should be noted that the formula for calculating the rate of change of elevation in the x-direction is: ,in, - This represents the elevation values ​​of the grid cells within a 3×3 window. The cells are arranged such that, centered on the cell whose slope is to be calculated, the nine cells within the 3×3 window are numbered in a row-column order from left to right and from top to bottom. Given the DTM raster size, the formula for calculating the rate of change of elevation in the y-direction is: .

[0119] The terrain relief is the difference between the maximum and minimum elevation values ​​of the grid cells within the window; the surface cutting depth is the difference between the average and minimum elevation values ​​of the grid cells within the window. The formula for calculating the elevation variation coefficient is: ,in, The standard deviation of the elevation of the grid within the window. This represents the average elevation value of the grid within the window. A value > 0.01 indicates uneven terrain. A value < 0.005 indicates that the terrain is flat.

[0120] To avoid missing micro-topographical scales with a single window, three windows with different pixel scales are selected to extract topographical factors. The specific settings are as follows:

[0121] The window sizes were selected as 31×31 pixels, 41×41 pixels, and 51×51 pixels, corresponding to actual field ranges of 1.55m×1.55m, 2.05m×2.05m, and 2.55m×2.55m, respectively. These values ​​were based on experimental results. The elevation variance converged within the 30-60 pixel window range. The 31-51 pixel window can cover the micro-topographic range of furrows and depressions within cultivated land, avoiding the omission of small-scale furrow factors by a single window.

[0122] For each farmland plot, five topographic factors are extracted under three windows to generate 15 sets of factor data.

[0123] The five initial topographic factors (elevation, slope, topographic relief, surface incision depth, and elevation variation coefficient) may exhibit multicollinearity (i.e., a high degree of linear correlation between the factors), leading to distorted parameter estimation in subsequent models. The variance inflation factor method is used to screen the multi-window topographic factors to obtain key topographic factors. The specific implementation is as follows:

[0124] The variance inflation factor (VIF) quantifies the degree of collinearity of a single factor due to the combined influence of all other factors, making it more suitable for multi-factor redundancy screening. The VIF is used to calculate the multiple correlation coefficient between a factor and all other factors. To measure collinearity, the formula for calculating the variance inflation factor is:

[0125]

[0126] in, Let be the variance inflation factor of the i-th factor. The regression determination coefficient is calculated when the i-th factor is the dependent variable and all other factors are independent variables.

[0127] It should be noted that the physical meaning of this formula is as follows: When = 1, the i-th factor is completely independent of the other factors. A value >10 indicates severe multicollinearity, and this factor should be removed. The larger the value, the stronger the linear influence of other factors on the factor, and the higher the redundancy.

[0128] For each factor, a regression model is constructed with all other factors as independent variables. The model is then fitted using the least squares method to obtain the regression results for each factor. Substitute into the formula to calculate Factors with VIF > 10 (surface cutting depth, elevation variation coefficient) are first removed because they are significantly affected by the linearity of other factors and have high information redundancy. The VIF values ​​of the remaining factors are recalculated (at this time, the number of factors is reduced, and it is necessary to verify whether collinearity still exists). Factors with VIF < 10 (the industry common threshold) are retained and recorded as key terrain factors. Finally, the key terrain factors are obtained: elevation, slope and terrain relief.

[0129] It should be noted that both the surface cutting depth and the topographic relief are based on "elevation difference". The elevation variation coefficient depends on elevation data. The two have "information overlap" with other factors, which inevitably leads to VIF>10. Elevation (benchmark), slope (local), and topographic relief (overall) describe the terrain from different dimensions, with no information redundancy, resulting in VIF<10. These three factors are obtained by quantitative screening based on real data.

[0130] S4: Based on the key terrain factors, the entropy weight method is used to perform weighted calculation to obtain the comprehensive membership degree of the multi-window terrain factors. The comprehensive membership degree of the multi-window terrain factors is compared with a preset threshold to select target slope points. The target slope points are spatially superimposed with cultivated land plots to obtain the preliminary field slope of the cultivated land plots.

[0131] Based on the aforementioned key terrain factors, an entropy weighting method is used for weighted calculation to obtain the comprehensive membership degree of the multi-window terrain factors. The specific implementation is as follows:

[0132] The three terrain factors extracted in step 3 have different dimensions. Directly using them in the calculation would lead to weight bias due to differences in their numerical ranges. Therefore, "linear normalization" is used to transform all factors into the [0,1] interval, laying the foundation for subsequent entropy weighting. The specific implementation is as follows:

[0133] The formula for calculating numerical normalization is:

[0134]

[0135] in, Let be the normalized value of the i-th terrain factor in the k-th pixel. Let be the original value of the i-th terrain factor in the k-th pixel. Let be the minimum value of the i-th topographic factor among all cultivated land pixels in the quadrat. It represents the maximum value of the i-th topographic factor among all cultivated land pixels in the quadrat.

[0136] After normalization, all factors are unified to the [0,1] interval, eliminating dimensional differences.

[0137] Traditional slope calculations often use equal weights or subjective weighting, which can easily lead to biased results. This method quantifies the information content of factors through their "information entropy." The lower the information entropy (the greater the factor's variability and the more information it contains), the greater its weight, thus achieving objective weighting. The specific implementation steps are as follows:

[0138] First, calculate the normalized percentage of each factor in each pixel, reflecting the contribution of that pixel's factor value to the overall picture. The formula for calculating the factor weight is:

[0139]

[0140] in, Let be the proportion of the i-th factor in the k-th pixel. Let be the normalized value of the i-th terrain factor in the k-th pixel. This represents the total number of pixels in all cultivated land areas within a sample plot.

[0141] Information entropy is a core indicator for measuring the amount of information in a factor. The formula for calculating information entropy is:

[0142]

[0143] in, Let i be the information entropy of the i-th factor. The standardized coefficient is... Let represent the proportion of the i-th factor in the k-th pixel.

[0144] The weights reflect the contribution of factors to the field slope, and the formula for calculating factor weights is as follows:

[0145]

[0146] in, Let i be the weight of the i-th factor. This represents the information utility value of the i-th factor; a larger value indicates that the factor contains more information. The sum of the information utility values ​​of all factors. The number of key terrain factors.

[0147] It should be noted that topographic relief has the highest weight, indicating that it has the greatest impact on field slope; slope has the lowest weight because the slope value of a locally flat area cannot reflect the overall micro-topography, which is consistent with the actual characteristics of farmland micro-topography.

[0148] Step 3 uses three windows (31×31, 41×41, and 51×51 pixels) to extract factors. It's necessary to fuse information from multiple windows to avoid omissions by a single window. Window weights are determined through linear weighting, and a comprehensive membership model is then constructed to screen target slope points (excluding interference from straw residues, fine soil clods, etc.). The specific implementation is as follows:

[0149] Window weights are determined based on the proportion of pixels in the window, and the formula for calculating window weights is:

[0150]

[0151] in, Let the weight of the j-th window be... Let J be the total number of pixels in the j-th window. This represents the total number of pixels in a 31x31 window. This represents the total number of pixels in a 41x41 window. This represents the total number of pixels in a 51x51 window.

[0152] For each window, a product-sum membership model is used. The formula for calculating the membership degree of a single window is as follows:

[0153]

[0154] in, Let k be the membership degree of the k-th pixel in each window. Let i be the weight of the i-th factor. It is the normalized value of the i-th factor in the k-th pixel.

[0155] The formula for calculating the comprehensive membership degree of multi-window terrain factors, which integrates the membership degrees of the three windows, is as follows:

[0156]

[0157] in, The multi-window terrain factor comprehensive membership degree of the k-th pixel. Let be the weight of the j-th window.

[0158] Grouping is determined based on the total number of pixels to ensure coverage of all membership distributions. The formula for calculating the number of hierarchical groups is as follows: ,in, The number of elevation extreme points within the quadrat is given, and then the class interval is calculated to reflect the interval of membership for each class: ,in, This represents the maximum value of the comprehensive membership degree of the multi-window terrain factors. To find the minimum value of the comprehensive membership degree of multi-window terrain factors, a membership degree frequency histogram was plotted. A "cliff-like drop" was observed at 0.627 (a sharp decrease in the number of interfering pixels, with the target pixel becoming dominant). Therefore, a threshold of 0.627 was determined. The comprehensive membership degree of multi-window terrain factors was compared with the preset threshold to filter out target slope points, and the comprehensive membership degree was retained. For pixels with a resolution greater than 0.627, remove interfering pixels such as stubble and small soil clods.

[0159] The original DTM elevation values ​​for depressions within the sample plots were too low due to the surrounding terrain obstruction, leading to an underestimation of the slope calculation. The slope of the depressions was corrected using a "reverse terrain flip" method, implemented as follows:

[0160] Using the average elevation of the sample plot DTM as a benchmark, an inverse terrain model is constructed by flipping the plot. The formula for calculating the elevation value of the raster in the inverse terrain model is as follows:

[0161]

[0162] in, These are the elevation values ​​of the raster in the inverse terrain model. This represents the average elevation value of the DTM quadrat. These are the elevation values ​​of the raster in the original DTM.

[0163] It should be noted that the significance of this formula is that depressions in the original DTM become "bumps" in the RDTM, eliminating the elevation underestimation caused by occlusion.

[0164] Repeat the above steps to extract three key factors (elevation, slope, and topographic relief) from the RDTM, perform normalization, entropy weighting, and multi-window comprehensive membership calculation, and screen the target slope points of the depressions: the slope value of the depression in the RDTM is equal to the actual slope value of the depression in the original DTM (because flipping only changes the elevation direction and does not change the slope magnitude), and directly extract the slope value screened in the RDTM as the corrected slope of the original depression.

[0165] The selected target slope points (including the depression slope corrected by RDTM) are spatially overlaid with the cultivated land plots to obtain the precise field surface slope for each plot. The specific implementation is as follows:

[0166] Using the "weighted average method" (weighted by pixel area to avoid edge pixel bias) to calculate slope, the preliminary formula for calculating the surface slope of a certain cultivated land plot is as follows:

[0167]

[0168] in, This represents the preliminary surface slope of a certain plot of farmland. Let be the slope value of the k-th target pixel within the plot. Let the area of ​​the k-th pixel be . This represents the total number of target pixels within the plot.

[0169] S5: Based on the coarse-resolution digital elevation model data, the area-weighted average method is used to calculate the topographic slope of the cultivated land plot; the K-means clustering method is used to calculate the main slope aspect correction coefficient of the cultivated land plot; the preliminary field slope, topographic slope, main slope aspect correction coefficient and shape index value are fitted by regression analysis using a power function model to obtain the field slope estimation model.

[0170] Based on the coarse-resolution digital elevation model data, the topographic slope of each cultivated land plot is extracted. The topographic slope of each plot is calculated using the slope formula in step 3 and the area-weighted average method in step 4. This ensures that the results reflect the overall topographic features of the site.

[0171] This step takes the preliminary field slope of the cultivated land plot generated in step 4 as the core, and combines the extracted cultivated land shape index (SI), the terrain slope of the coarse resolution digital elevation model and the newly added "main slope aspect correction coefficient" to construct an improved power function estimation model. By introducing terrain anisotropy correction (main slope aspect optimization), the problem of insufficient adaptability of traditional models in complex terrain areas is solved.

[0172] The model input parameters need to cover three dimensions: "topographic features, land parcel morphology, and slope aspect correction". Based on the data obtained from steps 2-4, standardization processing is performed to lay the foundation for model fitting. The specific implementation is as follows:

[0173] Choose 4 core parameters: dependent variable : Preliminary field slope of the sample plots of cultivated land obtained in step S4; Independent variable 1: Topographic slope Initial slope calculated based on a preprocessed 30m resolution SRTM DEM; Independent variable 2: shape index Correction factor for main slope aspect It is necessary to perform main slope aspect analysis and calculation to correct the influence of topographic anisotropy on slope.

[0174] To avoid model coefficient bias caused by differences in the magnitude of independent variables, Z-score standardization is used to convert all independent variables into standard values ​​with a mean of 0 and a standard deviation of 1. Z-score standardization can preserve the parameter distribution characteristics and is suitable for the preprocessing of independent variables in power function models.

[0175] Traditional models do not consider the differences in slope gradient across different slope aspects, leading to estimation errors in complex terrain areas. K-means clustering is used to calculate the main slope aspect correction coefficient for the digital terrain model, yielding the corrected coefficient for the main slope aspect of cultivated land plots. The specific implementation is as follows:

[0176] Extract the aspect value of each pixel from the high-precision DTM (5cm resolution) of the quadrat generated in step 2. The calculation formula is as follows:

[0177] in, For pixel slope value, The elevation change rate in the x-direction is... Let be the rate of change of elevation in the y-direction.

[0178] It should be noted that 0° is due north, 90° is due east, 180° is due south, and 270° is due west.

[0179] The K-means clustering algorithm is used to aggregate discrete aspect values ​​into "dominant aspect" (representing the dominant topographic orientation within the sample plot), and the angular aspect is converted to unit circle coordinates to avoid the continuity problem between 0° and 360°. The formula for calculating the aspect coordinate conversion is as follows:

[0180]

[0181] in, The slope angle, The coordinates are for the unit circle.

[0182] Based on the slope aspect distribution characteristics, eight initial centers are set (corresponding to eight directions: North, Northeast, East, Southeast, South, Southwest, West, and Northwest). The Euclidean distance from each slope aspect point to the eight centers is calculated, and the point is assigned to the nearest center. The center coordinates are updated, and this process is repeated until the center stabilizes. The number of iterations is 50, and the convergence threshold is 0.01. The formula for calculating the Euclidean distance is:

[0183]

[0184] in, and Let d be the coordinates of the d-th cluster center. and For the first The coordinates of the slope aspect point The Euclidean distance from each slope point to the 8 centers. This is the slope aspect point number.

[0185] Select the cluster center with the largest number of samples and convert its coordinates back to angular slope: , This reflects the angle between the long side of the cultivated land plot and the main slope direction. The smaller the angle (the smaller the plot extends along the main slope direction), the greater the correction required for slope estimation. (Main slope direction correction factor) The calculation formula is: ,in, The angle of the longer side of the plot. The direction of the main slope of the quadrat.

[0186] By performing regression analysis on the field slope, topographic slope, main slope aspect correction coefficient, and shape index value using a power function model, a field slope estimation model is obtained. The calculation formula for the field slope estimation model is as follows:

[0187]

[0188] in, This represents the preliminary surface slope of a certain plot of farmland. For terrain slope, The shape index, The correction factor for the main slope aspect. and These are the model coefficients, obtained through parameter fitting.

[0189] It should be noted that in step S4 It is calculated based on high-precision data within the sample plots, focusing on eliminating local interference, and is suitable for accurate estimation within a small area. Its purpose is to provide a reliable true value basis for fitting the field slope estimation model. (Step S5) It is the input variable of the field slope estimation model, and is used as the dependent variable along with macroscopic topographic parameters ( By fitting parameters such as slope, the slope that was accurately estimated in a small area can be extended to the entire area.

[0190] To facilitate the calculation of coefficients, the natural logarithm of both sides of the formula is taken, transforming it into a linear model, let... , , The model then simplifies to The least squares regression method was used to fit A and b. The 12 plot samples were divided into training set and validation set. To ensure sufficient samples for model training, the training set contained 8 samples and the validation set contained 4 samples, with a ratio of 2:1. This ensured that the training set covered different slope ranges. The coefficients A and b were calculated and substituted into the coefficients to obtain the field slope estimation model.

[0191] S6: Input the shape index value, topographic slope and main slope aspect correction coefficient of the area to be estimated into the field slope estimation model, and output the preliminary field slope estimation results.

[0192] The area to be estimated is large (e.g., 50 km²). 2 This step requires replacing the visual interpretation in step 3. An automated method is used to extract the boundaries of cultivated land plots to provide basic units for batch estimation. The Gaofen-7 satellite imagery obtained in step 1 is preprocessed, and water bodies in the image are selected as dark targets (reflectivity ≈ 0). Atmospheric aerosol optical thickness (AOD) is calculated and substituted into the 6S radiative transfer model to eliminate the influence of atmospheric scattering and absorption on the image reflectivity. The Gaofen-7 satellite imagery is then subjected to aerial triangulation. Based on the exterior orientation elements and SRTM DEM after aerial triangulation, orthorectification is performed using bilinear interpolation. The planar error of the corrected image is < 0.5m.

[0193] The Normalized Difference Vegetation Index (NDVI) was calculated using the near-infrared and red bands of Gaofen-7 satellite imagery to highlight the spectral differences between cultivated and non-cultivated land. The calculation formula is as follows:

[0194]

[0195] in, Normalized Difference Vegetation Index (NDVI) For near-infrared reflectivity, This refers to the reflectivity in the red band.

[0196] By statistically analyzing the NDVI of sample plots, a threshold for farmland identification was set at NDVI > 0.3. The farmland area was initially extracted, and the boundaries of the extracted farmland area were refined using the Canny operator. Areas with NDVI > 0.3 were set to 1 (farmland), and the rest were set to 0. A binary image was generated, and a 3×3 rectangular structuring element was used for closing operations to fill small holes (area < 10m²) inside the farmland. 2 To ensure the integrity of the land parcel, refer to the edge optimization technique in step 3, extract continuous farmland boundaries, convert edge pixels into vector lines, eliminate dangling nodes (length < 0.5m) through "topology check", and generate closed farmland parcel vector surfaces.

[0197] Based on the interpreted farmland plot vectors, the three input parameters required for the model are extracted in batches (topography slope, shape index, and main slope aspect correction coefficient). The parameter definitions are the same as in step 5, but need to be adapted to the characteristics of large-scale data. The specific implementation is as follows:

[0198] Using the slope formula from step 3, the overall slope is calculated in ArcGIS using a 30m×30m window to obtain slope raster data. Based on the vector surface of cultivated land plots, the slope raster is then divided into zones for statistical analysis. An area-weighted average is used to calculate the slope for each plot. This ensures that the results reflect the overall topographic features of the site.

[0199] For all arable land parcels, the arable land parcel shape index is calculated in batches using the formula in step 3. Referring to the main slope aspect analysis method in step 5, but expanding the scope to the entire area to be estimated, the global slope aspect is extracted from the 30m DEM. The clustering process in step 5 is followed, but the number of cluster centers is increased to 16 (to accommodate greater terrain diversity). Finally, the main slope aspect and main slope aspect correction coefficient of the area to be estimated are determined. .

[0200] The shape index value, topographic slope and main slope aspect correction coefficient of the area to be estimated are input into the field slope estimation model, and the preliminary field slope estimation results are output to realize the batch estimation of field slope.

[0201] S7: Acquire photon point cloud data of the area to be estimated, process the photon point cloud data through a filtering algorithm with grid continuity constraints to obtain the surface elevation control points of the area to be estimated, construct an error correction model based on the surface elevation control points, input the preliminary field slope estimation results into the error correction model, and output the farmland surface slope estimation results, thereby realizing the estimation of farmland surface slope.

[0202] Acquire ICESat-2 satellite ATL03-level photon point cloud data covering the sample plot area. Specific parameters and selection criteria are as follows:

[0203] Select ATL03 level (raw photon point cloud). This level of data contains raw photon coordinates, signal confidence levels, and other information, which can meet the requirements for subsequent filtering and extraction of surface elevation points. The distance error along the track should be <1m to ensure the spatial matching accuracy between the photon point cloud and the image data. Only photon points with a confidence level ≥3 are retained, while solar background noise points are removed. After initial screening, the photon point density within the sample plot area should be ≥1 point / m. 2 .

[0204] The photon point cloud within the sample plot area was divided into 15m×15m grids. The number of photon points in each grid was counted. Grids with less than 3 points were removed (judged as isolated noise points), and grids with 3 or more points were retained (judged as potential signal point clusters). After preliminary denoising, the photon point retention rate was >85%.

[0205] Large-scale estimations are affected by atmospheric correction errors in satellite imagery and DEM elevation deviations, resulting in systematic errors (overall overestimation or underestimation). Photon point cloud-assisted correction technology, utilizing ICESat-2 photon point cloud data, corrects these systematic errors. The specific implementation is as follows:

[0206] The photon point cloud data is processed using a filtering algorithm constrained by raster continuity to obtain the surface elevation control points of the area to be estimated. The initial raster size is calculated using the following formula:

[0207]

[0208] in, For grid width, This represents the highest elevation of the photon point within the sample plot area. This represents the lowest elevation of the photon point within the sample plot area.

[0209] Photon point clouds within the sample plot area are divided according to the initial grid size. The number of photon points in each grid is counted. Grids with ≥3 points are retained (determined as potential signal grids), and grids with <3 points are removed (determined as isolated noise grids). The connected components of all potential signal grids are labeled using the eight-connectivity criterion. The formula for calculating the grid continuity score is as follows:

[0210]

[0211] in, Raster continuity score; a higher score indicates stronger continuity between the raster and the surrounding terrain. The number of rasters in the connected region to the left of the raster to be detected. This represents the number of rasters within the connected region to the right of the raster to be detected.

[0212] Retain continuity score Grids with a value of ≥5 (defined as "signal grids") are discarded. For grids with a size <5 (identified as "noise grids"), an iterative scaling strategy is adopted to further improve filtering accuracy. In each iteration, the grid size is reduced to half of the previous size (scaling ratio r=2), and the grid division-continuity detection process is repeated. The iteration terminates when the grid size is ≤0.5m.

[0213] Fifty verification points (covering different slope grades) were evenly selected in the area to be estimated. An error correction model was constructed based on the surface elevation control points. The preliminary field slope estimation results were input into the error correction model, and the farmland field slope estimation results were output.

[0214] Estimated slope of the plot at each point The calculated slope of the field at 50 verification points, obtained through the preceding steps, was measured using RTK (accuracy ±0.1°) and used as the true value. Calculate the estimation error for each verification point. ,Think Independent variable, with With as the dependent variable, the formula for calculating the error correction model is:

[0215]

[0216] in, To estimate the error, Estimate the slope of the plot at each point. , and The correction coefficients are obtained through least squares regression fitting. =0.35, =-0.02, =0.005.

[0217] The error correction model is applied to the estimated slope of all plots, and the calculation formula is as follows:

[0218]

[0219] in, The results of the farmland surface slope estimation. To estimate the error.

[0220] This method acquires image data and a coarse-resolution digital elevation model of the area to be estimated. First, it refines the data to construct a high-precision digital ground model and orthophoto map of the sample plots. Then, it optimizes the farmland plots through visual interpretation, extracts topographic factors via a multi-window method (using the variance expansion factor method to remove redundancy), and uses entropy weighting to select target slope points. Finally, it combines an inverse topography model to correct the slope of depressions and obtain the true value of the field slope. Subsequently, it introduces a main slope aspect correction coefficient to construct an improved power function model. Finally, it automatically interprets farmland based on optical satellite imagery and corrects for errors in photon point clouds, achieving accurate and large-scale estimation of farmland slope.

[0221] This technical solution integrates multi-source data and constructs a high-precision digital ground model of sample plots. It generates dense point clouds through UAV motion recovery structure and multi-view stereo vision technology, and combines photon point cloud preliminary denoising. It can completely preserve micro-topographic features such as ridges and depressions in cultivated land, solving the problem that traditional single data sources are difficult to capture fine terrain and have insufficient accuracy, and providing reliable basic data support for subsequent slope calculation.

[0222] Based on the constructed high-precision digital ground model, the variance expansion factor method is used to quantify the multicollinearity of topographic factors in cultivated land plots to eliminate redundancy. Then, the entropy weight method is used to objectively calculate the weights of key topographic factors, avoiding the interference problem of topographic factors caused by traditional subjective weighting. This makes the selection and weight allocation of topographic factors more in line with the actual topographic characteristics of cultivated land, accurately distinguishes the influence of target slope points and interference features, and further improves the accuracy of field slope calculation within the sample plot.

[0223] K-means clustering is introduced to determine the main slope aspect of cultivated land and calculate the correction coefficient. An improved power function model is constructed, which breaks through the limitations of traditional models that only consider topographic slope and shape index and ignore topographic anisotropy. It can better adapt to the estimation of cultivated land slope in complex terrain areas, and make the model's slope prediction for different slope aspects more consistent with the actual tilt of the cultivated surface.

[0224] Based on the constructed power function model, the surface elevation control points are obtained by using the continuity constraint filtering of photon point cloud raster. A quadratic polynomial error correction model is constructed to effectively eliminate the accumulation of systematic errors in large-scale estimation. This solves the problem that the accuracy of the model is prone to decrease due to data deviation when applied to the whole domain, and further ensures the reliability of large-scale farmland slope estimation.

[0225] It should be noted that, in this document, relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, the phrase "comprising an element defined as..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0226] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A remote sensing-based method for estimating the slope of a ploughed field, characterized in that: Includes the following steps: Step S1: Obtain image data and coarse-resolution digital elevation model data for the area to be estimated; Step S2: Rasterize and edit the image data to obtain a digital ground model and orthophoto map; Step S3: Visually interpret the orthophoto map to obtain the distribution of cultivated land plots within the sample plots of the area to be estimated and calculate the shape index value of the cultivated land plots. Based on the digital terrain model, use the multi-window extraction method and slope formula to calculate and obtain multi-window terrain factors. Use the variance expansion factor method to screen the multi-window terrain factors and obtain key terrain factors. Step S4: Based on the key terrain factors, the entropy weight method is used to perform weighted calculation to obtain the comprehensive membership degree of the multi-window terrain factors. The comprehensive membership degree of the multi-window terrain factors is compared with a preset threshold to select target slope points. The target slope points are spatially superimposed with cultivated land plots to obtain the preliminary field slope of the cultivated land plots. Step S5: Based on the coarse-resolution digital elevation model data, the area-weighted average method is used to calculate the topographic slope of the cultivated land plots; the K-means clustering method is used to calculate the main slope aspect correction coefficient of the cultivated land plots using the digital ground model; and a power function model is used to perform regression analysis and fitting on the preliminary field slope, topographic slope, main slope aspect correction coefficient, and shape index value to obtain the field slope estimation model, including the following specific steps: The K-means clustering algorithm is used to aggregate discrete aspect values ​​into principal aspects, and the angular aspects are converted into unit circle coordinates to avoid the continuity problem between 0° and 360°. The calculation formula for aspect coordinate conversion is as follows: ； wherein is the slope angle, is the unit circle coordinate; Based on the slope aspect distribution characteristics, eight initial centers are set to correspond to eight directions: North, Northeast, East, Southeast, South, Southwest, West, and Northwest. The Euclidean distance from each slope aspect point to the eight centers is calculated, and the point is assigned to the nearest center. The center coordinates are updated, and this process is repeated until the center stabilizes. The number of iterations is 50, and the convergence threshold is 0.

01. The formula for calculating the Euclidean distance is: ； wherein, and are the coordinates of the d-th cluster center, and are the coordinates of the d-th slope point, are the Euclidean distances of each slope point to the 8 centers, is the slope point index; Select the cluster center with the largest number of samples and convert its coordinates back to angular slope: , This reflects the angle between the long side of the cultivated land plot and the main slope direction. The smaller the angle, the greater the correction required for slope estimation. (Main slope direction correction factor) The calculation formula is: ,in, The angle of the longer side of the plot. Main slope of the quadrat; Step S6: Input the shape index value, topographic slope and main slope aspect correction coefficient of the area to be estimated into the field slope estimation model, and output the preliminary field slope estimation results; Step S7: Obtain photon point cloud data of the area to be estimated, process the photon point cloud data using a raster continuity constraint filtering algorithm to obtain the surface elevation control points of the area to be estimated, construct an error correction model based on the surface elevation control points, input the preliminary field slope estimation results into the error correction model, and output the farmland surface slope estimation results, thereby realizing the estimation of farmland surface slope. The photon point cloud data is processed using a filtering algorithm constrained by raster continuity to obtain the surface elevation control points of the area to be estimated. The initial raster size is calculated using the following formula: ； in, The grid width, This represents the highest elevation of the photon point within the sample plot area. This represents the lowest elevation of the photon point within the sample plot area; Using the 8-connectivity standard, the connected components of all potential signal grids are labeled. The formula for calculating the grid continuity score is as follows: ； in, Raster continuity score; a higher score indicates stronger continuity between the raster and the surrounding terrain. The number of rasters in the connected region to the left of the raster to be detected. This represents the number of grid cells within the connected region to the right of the grid cell to be detected. Estimated slope of the plot at each point The calculated slope of the field at 50 verification points was obtained through the preceding steps and was then measured using RTK (Real-Time Kinematics) to obtain the true values. Calculate the estimation error for each verification point. ,Think Independent variable, with With as the dependent variable, the formula for calculating the error correction model is: ； in, To estimate the error, Estimate the slope of the plot at each point. , and The correction coefficients are obtained through least squares regression fitting. =0.35, =-0.02, =0.005; The error correction model is applied to the estimated slope of all plots, and the calculation formula is as follows: ； in, The results of the farmland surface slope estimation. To estimate the error.

2. The method for estimating farmland surface slope based on remote sensing according to claim 1, characterized in that: The process of rasterizing and editing the image data to obtain a digital terrain model and orthophoto map includes the following specific steps: The dense point cloud of the UAV was converted into a digital ground model with a resolution of 0.05m using ordinary kriging interpolation. The interpolation formula is as follows: ； in, The estimated elevation values ​​for the raster to be interpolated. The number of known points in the interpolation neighborhood. For the first Interpolation weights for known points For the first Elevation values ​​of known points; By combining UAV orthophotos, the elevation points of the digital ground model corresponding to the grids with gray values ​​<50 are removed; the grids with elevation values ​​>0.5m higher than the surrounding terrain are also removed. The minor undulations caused by median filtering are smoothed. The digital ground model is visually checked to ensure that 5cm-scale furrows and 20cm-scale depressions are retained without omission. Finally, the digital ground model and the corresponding orthophoto map are generated.

3. The method for estimating farmland surface slope based on remote sensing according to claim 2, characterized in that: The process of visually interpreting orthophoto maps to obtain the distribution of cultivated land parcels within the sample plots of the area to be estimated includes the following steps: Visual interpretation was performed based on orthophoto maps, and the following visual interpretation standards were set: cultivated land in black soil areas showed parallel ridge and furrow textures, cultivated land in plain areas showed uniform crop row textures, and the gray value range of cultivated land in orthophoto maps was 80-150. The boundaries of cultivated land plots were manually drawn in ArcGIS 10.8 software to generate initial cultivated land vector surfaces. A total of 12 cultivated land plots were extracted within the sample plots to ensure that there were no omissions or redundant drawings. The initial boundary is optimized by first performing histogram equalization on the orthophoto image to enhance the grayscale difference between cultivated and non-cultivated land. Edge detection is then performed using the Canny operator, with the following operator parameters set: Gaussian filter standard deviation. =1.0, low threshold T1=50, high threshold T2=150; The detected farmland edge vectors are superimposed with the visually interpreted initial boundaries, and then fused using the minimum distance method to finally generate the farmland plot boundary locations.

4. The method for estimating farmland slope based on remote sensing according to claim 3, characterized in that: The calculation of the shape index value of cultivated land plots includes the following steps: The formula for calculating the shape index of each arable land parcel is: ； in, This is a farmland plot shape index, reflecting the regularity of the plots. When =1, it is a square. The larger the plot of land, the more irregular its shape. The perimeter of the land parcel boundary. This refers to the area of ​​the land parcel.

5. The method for estimating farmland slope based on remote sensing according to claim 4, characterized in that: The process of calculating multi-window terrain factors based on the digital terrain model using a multi-window extraction method and slope formula includes the following steps: The slope of each grid cell is calculated using a 3×3 moving window. The calculation formula is as follows: ； in, This represents the slope value for a single grid cell, in degrees. Let x be the rate of change of elevation in the x-direction. The elevation change rate in the y-direction; The terrain relief is the difference between the maximum and minimum elevation values ​​of the grid cells within the window; the surface cutting depth is the difference between the average and minimum elevation values ​​of the grid cells within the window. The formula for calculating the coefficient of variation of elevation is: ,in, The standard deviation of the elevation of the grid within the window. This represents the average elevation value of the grid within the window; A value > 0.01 indicates uneven terrain. A value < 0.005 indicates flat terrain; To avoid missing micro-topographical scales with a single window, three windows with different pixel scales were selected to extract topographical factors: the window scales were 31×31 pixels, 41×41 pixels, and 51×51 pixels, corresponding to actual field ranges of 1.55m×1.55m, 2.05m×2.05m, and 2.55m×2.55m, respectively. The values ​​were based on experimental results, and the elevation variance converged within the 30-60 pixel window range. The 31-51 pixel window can cover the micro-topographical range of furrows and depressions within cultivated land, avoiding the omission of small-scale furrow factors by a single window. For each farmland plot, five topographic factors are extracted under three windows to generate 15 sets of factor data.

6. The method for estimating farmland surface slope based on remote sensing according to claim 5, characterized in that: The process of filtering multi-window terrain factors using the variance inflation factor method to obtain key terrain factors includes the following steps: The formula for calculating the variance inflation factor is: ； in, Let be the variance inflation factor of the i-th factor. denoted as the regression determination coefficient when the i-th factor is used as the dependent variable and all other factors are used as independent variables; For each factor, a regression model is constructed with all other factors as independent variables. The model is then fitted using the least squares method to obtain the regression results for each factor. Substitute into the formula to calculate Factors with VIF > 10 are first removed; the VIF values ​​of the remaining factors are recalculated, and the factors with VIF < 10 are retained as key terrain factors, thus obtaining the key terrain factors.

7. The method for estimating farmland surface slope based on remote sensing according to claim 6, characterized in that: The process of calculating the comprehensive membership degree of multi-window terrain factors using the entropy weighting method based on the key terrain factors includes the following specific steps: For each window, a product-sum membership model is used. The formula for calculating the membership degree of a single window is as follows: ； in, Let k be the membership degree of the k-th pixel in each window. Let i be the weight of the i-th factor. This is the normalized value of the i-th factor in the k-th pixel; The formula for calculating the comprehensive membership degree of multi-window terrain factors, which integrates the membership degrees of the three windows, is as follows: ； in, The multi-window terrain factor comprehensive membership degree of the k-th pixel. Let be the weight of the j-th window.

8. The method for estimating farmland surface slope based on remote sensing according to claim 7, characterized in that: The process of using a power function model to perform regression analysis and fitting on the preliminary field slope, topographic slope, main slope aspect correction coefficient, and shape index values ​​to obtain a field slope estimation model includes the following specific steps: By performing regression analysis on the field slope, topographic slope, main slope aspect correction coefficient, and shape index value using a power function model, a field slope estimation model is obtained. The calculation formula for the field slope estimation model is as follows: ； in, This represents the preliminary surface slope of a certain plot of farmland. For terrain slope, The shape index, The correction factor for the main slope aspect. and These are the model coefficients, obtained through parameter fitting.

Images