Land area measurement method and system for land planning and design

Abstract

The invention discloses a land area measurement method and system for land planning and design, and the method comprises the steps: collecting a sequence aerial image of a target land parcel, and synchronously recording the positioning and attitude determination data of each frame of image; edge sharpening processing and multi-view three-dimensional reconstruction are carried out on the sequence aerial image, and a digital elevation model with texture features is generated; performing terrain gradient analysis and edge contour extraction based on a digital elevation model, and identifying a land block boundary through an adaptive threshold segmentation algorithm and generating an initial vector boundary diagram; performing topological relation verification and fragment pattern spot fusion on the initial vector boundary diagram, and outputting an optimized plot vector boundary; and carrying out coordinate registration on the optimized plot vector boundary and a geographic information system base map, calculating a projection area through integral operation, and generating a final measurement report containing a topographic relief coefficient. According to the embodiment of the invention, the operation efficiency, boundary recognition precision and result reliability of land area measurement under complex terrains can be improved.

Description

Technical Field

[0001] This invention belongs to the field of land surveying technology, and in particular to a method and system for measuring land area for land planning and design. Background Technology

[0002] In the field of land planning and design, accurate and efficient land area measurement is fundamental for land resource management, project planning, and ownership delineation. Traditional measurement methods mainly rely on manual ground mapping using total stations and RTK, which are labor-intensive, inefficient, and pose safety risks and blind spots in complex terrain or large areas. In recent years, UAV remote sensing technology has been introduced into this field, using aerial imagery for measurement. However, existing technologies are mostly limited to generating orthophotos for two-dimensional planar measurement, failing to fully consider the significant impact of actual terrain undulations on the projected area. When dealing with ambiguous plot boundaries and complex surface textures, these methods have limited automation and boundary recognition accuracy, are susceptible to interference from shadows and vegetation occlusion, leading to deviations in the generated vector boundaries, and the measurement results fail to effectively integrate three-dimensional terrain information, making it difficult to meet the actual needs of high-precision land planning and design for "surface area" rather than "projected area." Summary of the Invention

[0003] The purpose of this invention is to provide a land area measurement method and system for land planning and design, so as to overcome the shortcomings of the prior art and improve the operational efficiency, boundary identification accuracy and reliability of land area measurement in complex terrain.

[0004] One embodiment of this application provides a method for measuring land area for land planning and design, the method comprising: A sequence of aerial images of the target site is collected by a drone equipped with a multispectral sensor, and the positioning and attitude data of each frame of the image are recorded simultaneously. The sequence of aerial images is subjected to edge sharpening and multi-view 3D reconstruction to generate a digital elevation model with texture features. Based on the digital elevation model, terrain gradient analysis and edge contour extraction are performed, and the plot boundaries are identified and an initial vector boundary map is generated by an adaptive threshold segmentation algorithm. The initial vector boundary map is subjected to topological relationship verification and fragment patch fusion. The elevation error is corrected by a dynamic weight compensation algorithm, and the optimized land parcel vector boundary is output. The optimized land parcel vector boundary is registered with the geographic information system base map, the projected area is calculated through integration, and a final measurement report including the terrain undulation coefficient is generated.

[0005] Another embodiment of this application provides a land area measurement system for land planning and design, the system comprising: The acquisition module is used to acquire a sequence of aerial images of the target site using a drone equipped with a multispectral sensor, and to simultaneously record the positioning and attitude data of each frame of the image. The processing module is used to perform edge sharpening and multi-view 3D reconstruction on the sequence of aerial images to generate a digital elevation model with texture features. The analysis module is used to perform terrain gradient analysis and edge contour extraction based on the digital elevation model, identify land parcel boundaries and generate an initial vector boundary map through an adaptive threshold segmentation algorithm; The verification module is used to perform topological relationship verification and fragment patch fusion on the initial vector boundary map, correct elevation errors using a dynamic weight compensation algorithm, and output the optimized land parcel vector boundary. The registration module is used to register the optimized land parcel vector boundary with the geographic information system base map, calculate the projected area through integration, and generate a final measurement report including the terrain undulation coefficient.

[0006] Another embodiment of this application provides a storage medium storing a computer program, wherein the computer program is configured to execute the method described in any of the preceding claims when running.

[0007] Another embodiment of this application provides an electronic device including a memory and a processor, wherein the memory stores a computer program and the processor is configured to run the computer program to perform the method described in any of the preceding claims.

[0008] Compared with existing technologies, the land area measurement method for land planning and design provided by this invention can improve the operational efficiency, boundary identification accuracy and reliability of land area measurement in complex terrain. Attached Figure Description

[0009] Figure 1 A hardware structure block diagram of a computer terminal for a land area measurement method for land planning and design provided in an embodiment of the present invention; Figure 2 A flowchart illustrating a land area measurement method for land planning and design provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of a land area measurement system for land planning and design, provided as an embodiment of the present invention. Detailed Implementation

[0010] The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0011] The present invention first provides a method for measuring land area for land planning and design. This method can be applied to electronic devices, such as computer terminals, specifically ordinary computers.

[0012] The following detailed explanation uses a computer terminal as an example. Figure 1 This is a hardware structure block diagram of a computer terminal for a land area measurement method for land planning and design, provided as an embodiment of the present invention. (See diagram below.) Figure 1 As shown, the computer device includes a processor, memory, and network interface connected via a system bus, wherein the memory may include non-volatile storage media and internal memory.

[0013] See Figure 2 The present invention provides a method for measuring land area for land planning and design, which may include the following steps: S201 uses a drone equipped with a multispectral sensor to collect a sequence of aerial images of the target site and simultaneously records the positioning and attitude data of each frame of the image. Specifically, the flight path of the drone can be planned according to the scope and shape of the target site, and the flight path plan and sensor parameter configuration can be determined; First, the boundary vector data of the target plot needs to be imported through Geographic Information System (GIS) software (such as ArcGIS) to clarify the spatial range and geometric shape of the plot. Taking a rectangular target plot as an example, its latitude and longitude range is 118.5000°-118.5200° east longitude and 30.2000°-30.2150° north latitude, with an area of ​​about 50 acres. The long side extends along the east-west direction, and the short side extends along the north-south direction, with no obvious irregular protrusions.

[0014] Flight path planning uses drone ground station software (such as DJI GS Pro), selecting the "grid flight" mode (suitable for full coverage data collection of regular plots). Path parameters need to be calculated in conjunction with sensor swath width and overlap rate. Flight altitude: Set according to the ground resolution requirements of the multispectral sensor. If a sampling accuracy of 5cm / pixel is required (to meet the error requirement of ≤0.5% for land area measurement), and the sensor selected is MicaSenseAltum (focal length 12mm, pixel size 3.75μm), then the flight altitude H = (ground resolution × focal length) / pixel size = (0.05m × 12 × 10) / 2000m. -3 m) / (3.75×10 - 6 m)=160m, but we actually used 150m to reserve redundancy; Line spacing: The sensor's horizontal field of view is 60°, and the ground swath width W = 2 × H × tan(field of view / 2) = 2 × 150 × tan(30°) ≈ 173.2m. To ensure an 80% lateral overlap rate between adjacent line images (to avoid missing shots), the line spacing D = W × (1 - overlap rate) = 173.2 × 0.2 ≈ 34.6m, rounded to 35m. Flight speed: The acquisition frequency of the multispectral sensor is set to 2Hz (1 frame is acquired every 0.5 seconds). To ensure a longitudinal overlap rate of 70%, the flight speed V = ground resolution × acquisition frequency / (1 - longitudinal overlap rate) = 0.05 × 2 / 0.3 ≈ 0.33 m / s. The actual speed was adjusted to 5 m / s (balancing efficiency and stability, the acquisition frequency was increased to 10Hz for adaptation, ensuring that the overlap rate meets the standard). Flight direction: Fly along the long side of the plot (east-west direction) to reduce the impact of heading deviation on image stitching. The starting point is set at the northwest corner of the plot (118.5000°E, 30.2150°N) and the ending point is set at the southeast corner (118.5200°E, 30.2000°N). A total of 12 horizontal flight paths are planned to cover the entire plot.

[0015] Sensor parameter configurations must match land planning and surveying requirements: Band selection: The multispectral sensor selects four core bands: blue (450-510nm, to distinguish water bodies from vegetation), green (560-620nm, to reflect vegetation cover), red (660-730nm, to identify plot boundaries), and near-infrared (770-890nm, to distinguish soil types), and eliminates redundant bands to reduce data volume. Exposure parameters: Use automatic exposure mode, ISO fixed at 100 (to reduce noise), shutter speed range 1 / 1000-1 / 2000s (to avoid blurring caused by flight shake); Data format: Images are saved in TIFF format (preserving original spectral information, uncompressed), with each frame image being approximately 25MB in size. Positioning and orientation data are simultaneously saved in CSV format for easy subsequent processing.

[0016] Control the drone to fly according to the flight path plan, collect a sequence of aerial images through multispectral sensors, and record high-precision positioning and attitude data of each frame of the image to generate the original aerial images and positioning and attitude dataset; The DJI M300RTK drone was selected as the flight platform. This model supports centimeter-level positioning (equipped with a D-RTK2 high-precision positioning module, with horizontal positioning accuracy of ±1cm+1ppm and vertical positioning accuracy of ±2cm+1ppm) and can stably carry multispectral sensors for long-term flight.

[0017] Equipment debugging must be completed before flight: Positioning calibration: Set up a ground control point near the plot (known latitude and longitude 118.5100°E, 30.2075°N, elevation 50.2m), and perform differential positioning calibration with the control point using the UAV RTK module to ensure that the positioning data deviation is <2cm; Sensor warm-up: After the multispectral sensor is powered on, it should be warmed up for 10 minutes until the spectral response stabilizes (based on ground testing, the spectral reflectance fluctuation of the same target is <1%) before data acquisition begins. Path loading: Import the flight path planned by the ground station into the UAV flight control system, confirm that the parameters such as the inflection point, altitude, and speed of the flight path are correct, and set up safety mechanisms such as "return to home when contact is lost" and "return to home when low battery" (triggered when the remaining battery is 20%).

[0018] Data acquisition process during flight: The UAV takes off automatically from the takeoff point (open area on the west side of the plot), climbs to a flight altitude of 150m, and flies along the preset grid route. The multispectral sensor acquires a sequence of aerial images at a frequency of 10Hz. Each frame of the image records the acquisition timestamp through EXIF ​​information (accurate to milliseconds, such as 20251001103000123). The D-RTK2 module records positioning and attitude data at a frequency of 20Hz, including the latitude and longitude (e.g., 118.5105°E, 30.2080°N), geodetic elevation (150.3m, based on the WGS84 coordinate system), and attitude angles (roll angle 0.2°, pitch angle -0.1°, yaw angle 180.4°, used for viewpoint correction in subsequent 3D reconstruction) at the time of each image acquisition. The flight control system uses a hardware timestamp synchronization mechanism to bind each frame of image with the closest positioning and attitude data (time difference ≤ 20ms) to avoid data misalignment. For example, the positioning and attitude data corresponding to a certain frame of image (timestamp 20251001103000123) is "longitude 118.510523°, latitude 30.208045°, elevation 150.32m, roll 0.2°, pitch -0.1°, yaw 180.4°".

[0019] After the flight, a set of original aerial images (a total of 800 frames were collected, covering the plot and 10% of the surrounding area to avoid missing boundary information) and a positioning and attitude dataset (800 records, one-to-one correspondence with the images) were generated. The data was stored on the drone's built-in SD card (capacity 128GB, remaining space ≥50GB). There was no obvious overexposure or underexposure in the single frame images, and the positioning data was continuous and uninterrupted.

[0020] The original aerial images and positioning and attitude determination datasets are time-synchronized and quality-checked, blurry or missing data frames are removed, and the verified sequence of aerial images and synchronized positioning and attitude determination data are generated.

[0021] Time synchronization verification aims to ensure the temporal consistency between each frame of image and the corresponding positioning and pose data, avoiding deviations in subsequent 3D reconstruction due to hardware latency. Extract the acquisition timestamp (T_img) from the EXIF ​​information of the original aerial image and the timestamp (T_pos) from the positioning and attitude data, and calculate the time difference ΔT=|T_img-T_pos|; The synchronization threshold is set to 50ms (based on UAV hardware latency testing, the maximum latency between the sensor and the positioning module is 35ms). If ΔT > 50ms, it is considered out of sync, and the frame image and corresponding positioning data are discarded. If ΔT ≤ 50ms, synchronization is considered valid. For example, for a frame image T_img=20251001103000123, corresponding to T_pos=20251001103000160, ΔT=37ms≤50ms, so synchronization is valid; for another frame ΔT=62ms, synchronization fails, and it is discarded.

[0022] Quality checks are divided into image sharpness checks and data integrity checks to ensure that the data source meets the requirements for subsequent processing. Image sharpness check: The BRISQUE (Blind / Referenceless Image Spatial Quality Evaluator) algorithm is used to calculate the sharpness score of each frame (range 0-100, the lower the score, the higher the sharpness); a sharpness threshold of 30 is set (after testing, images with a score ≤30 can clearly extract the boundaries of ground features, while images with a score >30 have obvious blur). If a frame has a score of 35 (due to blurred edges caused by flight shaking), it is discarded. Data integrity check: Traverse the positioning and attitude determination dataset. If there is no positioning record within a certain period of time (such as 1 consecutive second) (such as GPS signal obstruction), it is determined to be data missing, and all images within the corresponding time period are removed. For example, if a certain flight route is obstructed by trees, the positioning data is missing for 0.8 seconds, and the corresponding 4 frames of images (acquisition frequency 10Hz) are all removed.

[0023] After verification, 12 frames were removed from the original 800 frames due to time synchronization issues, 23 frames were removed due to insufficient sharpness, and 5 frames were removed due to missing data. This resulted in a final sequence of 760 verified aerial images and 760 synchronized positioning and attitude data points, achieving a data validity rate of 95%. A verification report was also generated, recording the reasons for removals, the number of removals of each type, and statistical information on the remaining data (e.g., average sharpness score of 22, average time difference of 28ms), ensuring the reliability of the data source for subsequent 3D reconstruction.

[0024] S202, perform edge sharpening processing and multi-view 3D reconstruction on the sequence of aerial images to generate a digital elevation model with texture features; Specifically, radiometric correction and noise filtering preprocessing can be performed on the verified sequence of aerial images to eliminate illumination differences and sensor noise, and generate a preprocessed image sequence. Although blurred frames have been removed from the verified sequence of aerial images, two types of interference still exist: one is the radiation error of the multispectral sensor (such as dark current and gray value shift caused by atmospheric scattering), and the other is sensor readout noise (such as sporadic bright spots and gray value fluctuations). These need to be eliminated through radiation correction and noise filtering to ensure the accuracy of subsequent edge sharpening and 3D reconstruction.

[0025] 1. Radiation Correction Radiometric correction addresses grayscale distortion in multispectral images and is performed in two steps: Dark current correction: Multispectral sensors still generate dark current (inherent pixel noise) even in the absence of light. It is necessary to acquire "dark frames" (black field images with the lens covered) and calculate the average dark current. For example, the average dark current of a blue band dark frame from a certain MicaSenseAltum sensor, after multiple acquisitions, is 15DN (DN is a digital quantization value, ranging from 0-255). The DN value of each pixel in the original image is subtracted from the corresponding band's dark current value. For instance, an original blue band pixel of 200DN will be corrected to 185DN, eliminating the high grayscale caused by dark current. Atmospheric correction: A simplified atmospheric radiative transfer model (6S model) was adopted. The meteorological parameters at the time of shooting (atmospheric visibility 10km, from the on-site meteorological station; solar zenith angle 30°, calculated based on the shooting time of 10:30) and sensor parameters (flight altitude 150m, focal length 12mm) were input to correct the influence of atmospheric scattering on the blue band. The original blue band image had an overall gray value of 180DN due to atmospheric scattering. After correction, it was reduced to 150DN, which was consistent with the gray range of the green band (145DN) and the red band (155DN), eliminating the illumination difference of "blue band being too bright".

[0026] 2. Noise filtering Gaussian filtering was used to remove sensor readout noise (such as random fluctuations in grayscale values ​​of ±5 DN). A 3×3 filter kernel was selected (a kernel size that is too small cannot completely filter noise, while a kernel size that is too large will blur the boundaries of objects), with a standard deviation σ=1.0 (σ controls the filtering intensity; 1.0 can smooth noise while preserving edge details). Taking a certain red band image as an example, sporadic bright spots (grayscale value 220 DN, surrounding pixels 180 DN) caused by noise in the original image were filtered by Gaussian. The grayscale value of the bright spots became 185 DN, and the transition with the surrounding pixels was natural. At the same time, the grayscale difference (180 DN vs 150 DN) of the field ridge boundary remained at 30 DN and was not blurred.

[0027] The final preprocessed image sequence has a gray value fluctuation of ≤3% in each band (gray value difference of the same ground object in different frames ≤5DN), and the noise standard deviation is reduced from 8DN to 3DN, laying the foundation for subsequent edge enhancement.

[0028] The Laplacian edge sharpening algorithm is applied to the preprocessed image sequence to perform convolution operation, thereby enhancing the boundary features of ground objects and generating an edge-sharpened image sequence. The grayscale contrast of the preprocessed image features (such as field ridges, roads, and plot edges) is low (usually only 20-30 DN), so the Laplacian edge sharpening algorithm is needed to enhance the boundary details. The Laplacian operator is a second-order differential operator that can highlight areas with drastic grayscale changes (i.e., boundaries), and it has a significant enhancement effect on linear features (such as field ridges).

[0029] 1. Laplace operator selection A 3×3 Laplacian operator template is used: [[0,1,0],[1,-4,1],[0,1,0]]. This template has the strongest edge response in the cross direction (which matches the common linear feature boundary direction in land planning) and avoids noise amplification in the diagonal direction.

[0030] 2. Convolution operation process For each pixel (except edge pixels) of the preprocessed image, perform convolution as follows: Take the 3×3 neighboring pixels of pixel (x,y), and denote them as N(x-1,y-1), N(x-1,y), N(x-1,y+1), N(x,y-1), N(x,y), N(x,y+1), N(x+1,y-1), N(x+1,y), N(x+1,y+1); Multiply the neighboring pixels by the corresponding positions of the operator template and sum them to obtain the edge enhancement value: Enhancement value = N(x-1,y)×1+N(x,y-1)×1+N(x,y)×(-4)+N(x,y+1)×1+N(x+1,y)×1; To avoid overexposure at the edges, the enhancement value is proportionally added to the original pixel value: Sharpened pixel value = Original pixel value + 0.5 × Enhancement value (0.5 is the addition coefficient to balance the enhancement effect and the naturalness of the image).

[0031] 3. Example of enhanced effect Taking the field ridge boundary in the green band image as an example: After preprocessing, the grayscale difference between the field ridge pixel (160DN) and the adjacent farmland pixel (130DN) is 30DN; after convolution operation, the edge enhancement value = 130×1 + 130×1 + 160×(-4) + 130×1 + 130×1 = -120; after sharpening, the field ridge pixel value = 160 + 0.5×(-120) = 100DN, the farmland pixel value = 130 + 0.5×120 = 190DN, the grayscale difference is increased to 90DN, and the field ridge boundary changes from "blurry thin line" to "clear dark line".

[0032] The generated edge-sharpened image sequence improves the contrast of ground feature boundaries by 2-3 times without significant noise amplification (the absolute value of edge enhancement is controlled within 100 to avoid pixel overflow in the 0-255 range), providing clear boundary features for subsequent key point matching.

[0033] Based on synchronous positioning and attitude determination data, a scale-invariant feature transformation algorithm is used to match key points in multi-view images and generate a set of matching feature points. For multi-view aerial images (with a difference of 5-10° in shooting angle between adjacent frames and a scale difference of <2%), it is necessary to establish spatial relationships between images through stable key point matching. The Scale Invariant Feature Transform (SIFT) algorithm has scale invariance and rotation invariance, and can achieve reliable matching under complex viewpoint changes. Synchronous positioning and attitude data (such as camera position and attitude angle) can help filter out incorrect matches.

[0034] 1. Core Process of SIFT Algorithm Scale-space extremum detection: A Gaussian pyramid (6 groups, 4 layers per group) is constructed, with a scale factor of 1.2 between adjacent groups (controlling the scale variation step size). Gaussian difference images are obtained by subtracting adjacent layer images, and local extrema (potential keypoints) are detected in the difference images. For example, in a frame of red band image, 200 potential keypoints were detected in the difference images at scales σ=2.0 and σ=2.4. Key point localization: Secondary Taylor expansion is performed on potential key points to remove key points with low contrast (contrast threshold 0.04, to avoid noise points) and edge response (edge ​​threshold 10, to avoid unstable edge points), and finally 150 stable key points are retained. Their scale (e.g. σ=2.2) and image coordinates (e.g. (320,240)) are recorded. Direction assignment: Centered on the key point, calculate the direction and magnitude of the pixel gradient in a 16×16 neighborhood, generate a 36-bin (10° per bin) orientation histogram, and take the peak direction (e.g., 60°) as the main direction; if the difference between the secondary peak and the main peak is <20%, add an auxiliary direction (e.g., 150°) to enhance the matching robustness; Feature descriptor generation: The neighborhood of the key point is divided into 4×4 sub-blocks. The gradient histogram of each sub-block is calculated in 8 directions to generate a 128-dimensional feature descriptor (4×4×8=128). The descriptor elements are normalized to 0-1 to eliminate the influence of illumination (e.g., a key point descriptor fragment is [0.12,0.08,0.15,...,0.05]).

[0035] 2. Key point matching and filtering The matching uses the "nearest neighbor distance ratio method": for each key point in image A, find the nearest neighbor (smallest distance) and the second nearest neighbor (second smallest distance) key points in image B, and calculate the distance ratio (nearest neighbor distance / second nearest neighbor distance). If the ratio is <0.8 (experiments have verified that this threshold can control the false matching rate within 2%), it is considered a correct match.

[0036] Further filtering is performed using synchronous positioning and attitude data: the 3D distance between theoretical matching points is calculated based on camera extrinsic parameters (camera position and attitude angle of the two frames). If the 3D distance between the actual matching points (calculated by triangulation) deviates from the theoretical value by more than 0.5m, the matching pair is discarded. For example, in the SIFT matching of image A (frame 100, camera position 118.5100°E, 30.2075°N, 150m) and image B (frame 101, camera position 118.5105°E, 30.2076°N, 150m), 150 pairs of matching points are initially matched. After filtering by distance ratio and extrinsic parameters, 120 pairs of correct matching points are retained. The generated matching feature point set contains the coordinates (e.g., A image (320, 240) corresponds to B image (350, 250)), scale, and orientation information of each pair of matching points, providing sparse correspondence for 3D reconstruction.

[0037] A multi-view stereo vision algorithm is used to reconstruct a dense 3D point cloud from the matching feature point set, generating a preliminary 3D point cloud model. Sparse matching points (e.g., 120 pairs / frame) cannot reflect terrain details. Dense point clouds need to be generated by multi-view stereo vision algorithms (PMVS, Patch-Based Multi-View Stereo). The PMVS algorithm can expand from sparse points to dense patches, which is suitable for the measurement needs of details such as small slopes and field ridges in land planning (accuracy ≤ 0.1m).

[0038] 1. Core Process of PMVS Algorithm Patch initialization: Generate 3×3 square patches centered on SIFT matching points (patch size adapted to 1920×1080 image resolution, balancing detail and efficiency). Each patch determines its corresponding position in the multi-view image through matching relationships. Patch matching and verification: For each patch, find the patch with the highest grayscale correlation in other viewpoint images (using normalized cross-correlation coefficient NCC, threshold 0.7, the closer the NCC is to 1, the more reliable the match); combined with synchronous positioning and attitude data (camera intrinsic parameters: focal length 12mm, pixel size 3.75μm; extrinsic parameters: position and attitude angle of the camera in each frame), calculate the three-dimensional coordinates of the patch through triangulation - for example, the intersection of the projection rays of the patch on camera A (position 118.5100°E, 30.2075°N, 150m) and camera B (position 118.5105°E, 30.2076°N, 150m) is the three-dimensional coordinate (118.5102°E, 30.2075°N, 50.3m); Dense point cloud generation: For verification patches with NCC≥0.7, the three-dimensional coordinates of their centers are extracted as point cloud points; outliers are removed by statistical filtering (the distance between each point and its neighboring points is calculated, and points with a distance mean ±2 standard deviations are removed). The point cloud density is set to 10 points per square meter (to meet the accuracy of topographic surveying, approximately 330,000 points are generated for a 50-acre plot).

[0039] 2. Preliminary Point Cloud Model Example The preliminary 3D point cloud model coordinates of the target plot are: 118.5000°-118.5200° E, 30.2000°-30.2150° N, and elevation 48.5m-52.3m. The point cloud has no obvious holes (hole area <0.1%), and the point cloud is continuous in areas with a slight slope (slope 2°), clearly showing topographic details such as field ridges (elevation difference 0.3m) and low-lying areas (elevation 48.5m). Verified by ground control points (known elevation 50.2m), the point cloud elevation error is ≤0.05m, which meets the accuracy requirements of land surveying.

[0040] The texture information of the edge-sharpened image sequence is mapped onto the preliminary 3D point cloud model, and a digital elevation model with texture features is generated through surface reconstruction algorithms.

[0041] The initial point cloud model only contains three-dimensional coordinates. It is necessary to map the texture information (such as the color of ground features and boundary details) of the edge-sharpened image, and then generate a continuous digital elevation model (DEM) through surface reconstruction. The DEM is the core terrain data for land area measurement.

[0042] 1. Texture mapping Texture image selection: For each point cloud point, calculate the angle between the camera optical axis and the point cloud point's normal vector (the normal vector reflects the terrain's tilt direction) based on its 3D coordinates and camera extrinsic parameters. Select images with an angle <30° as texture sources (to avoid texture stretching). For example, if a point cloud point (elevation 50.3m) has the best viewing angle in a sharpened image at frame 100, select the texture from that image. Texture coordinate calculation: The 3D coordinates of the point cloud points are projected onto the texture image using camera intrinsic parameters to obtain pixel coordinates (texture coordinates). For example, the pixel coordinates (325, 245) of the point cloud point (118.5102°E, 30.2075°N, 50.3m) projected onto the image frame 100 are extracted. The multispectral texture (red band 150DN, green band 160DN, blue band 140DN) of this pixel and its 3×3 neighborhood is extracted. Texture attachment: The texture information is assigned to the point cloud points. The textures of adjacent point cloud points are blended through mean fusion. For example, the texture of the field ridge area smoothly transitions from dark texture (100DN) to bright texture of farmland (190DN) without obvious splicing marks.

[0043] 2. Surface Reconstruction and DEM Generation A Poisson surface reconstruction algorithm is used to generate a continuous terrain surface. This algorithm can fit a smooth surface using point cloud normal vectors. The parameter settings are as follows: Normal vector estimation radius: 0.5m (set according to point cloud density to ensure that the normal vector accurately reflects the terrain tilt); Reconstruction depth: 8 (controls the level of detail in the DEM; a depth of 8 can preserve 0.1m of terrain undulation). DEM resolution: 1m×1m (each grid cell represents a 1m×1m area on the ground, storing the average elevation).

[0044] The final generated digital elevation model contains 333×150 grid cells (covering a 50-acre plot), with an elevation range of 48.5m-52.3m. The texture perfectly matches the terrain—the texture of the field ridge area is linear dark texture (corresponding to the boundary in the sharpened image), and the texture of the farmland area is uniform bright texture. The elevation accuracy of the DEM has been verified, with an error of ≤0.05m from the ground control points, and it can be directly used for subsequent terrain gradient analysis and boundary extraction.

[0045] S203, Based on the digital elevation model, perform terrain gradient analysis and edge contour extraction, identify land parcel boundaries and generate an initial vector boundary map through an adaptive threshold segmentation algorithm; Specifically, terrain gradient calculations can be performed on digital elevation models, and the Sobel operator can be used to solve for the spatial derivative of elevation values ​​to generate terrain gradient maps. The Digital Elevation Model (DEM) is a 1m×1m grid of data generated in previous steps. Each grid cell corresponds to a 1m×1m area on the ground, storing the average elevation value (in meters) of that area. For example, the DEM grid elevation range for a target plot is 48.5m-52.3m. The grid elevation in the field ridge area is 0.2-0.3m higher than that of adjacent farmland, while the elevation in low-lying areas is concentrated between 48.5-49.0m. The terrain gradient reflects the rate of change of elevation values ​​in space. Areas with a high rate of change (such as field ridges and plot edges) are often the core locations of plot boundaries. Therefore, the gradient needs to be calculated using the Sobel operator to highlight these areas.

[0046] The Sobel operator is a commonly used edge detection operator. It calculates the spatial derivatives of the image in two orthogonal directions, x (horizontal, east-west) and y (vertical, north-south), to obtain the magnitude and direction of the gradient. The larger the gradient magnitude, the more drastic the elevation change at that location. In its implementation, the Sobel operator first defines the x-direction template (Sobel_x) and the y-direction template (Sobel_y). The Sobel_x template is [[-1,0,1],[-2,0,2],[-1,0,1]], used to detect elevation changes in the vertical direction (north-south); the Sobel_y template is [[-1,-2,-1],[0,0,0],[1,2,1]], used to detect elevation changes in the horizontal direction (east-west).

[0047] The calculation process centers on each grid cell (i,j) of the DEM and takes the elevation values ​​of 9 grid cells within its 3×3 neighborhood (denoted as Z(i-1,j-1), Z(i-1,j), Z(i-1,j+1), Z(i,j-1), Z(i,j), Z(i,j+1), Z(i+1,j-1), Z(i+1,j+1), Z(i+1,j+1)). These values ​​are then convolved with the Sobel_x and Sobel_y templates to obtain the gradients Gx and Gy in the x and y directions, respectively. Gx=Z(i-1,j-1)×(-1)+Z(i-1,j)×0+Z(i-1,j+1)×1+Z(i,j-1)×(-2)+Z(i,j)×0+Z(i,j+1)×2+Z(i+1,j-1)×(-1)+Z(i+1,j)×0+Z(i+1,j+1)×1; Gy=Z(i-1,j-1)×(-1)+Z(i-1,j)×(-2)+Z(i-1,j+1)×(-1)+Z(i,j-1)×0+Z(i,j)×0+Z(i,j+1)×0+Z(i+1,j-1)×1+Z(i+1,j)×2+Z(i+1,j+1)×1.

[0048] The gradient magnitude G is then calculated as √(Gx² + Gy²), where the unit of magnitude G is meters per meter (i.e., dimensionless, reflecting elevation changes per unit distance). For example, the elevation values ​​(in meters) of a 3×3 neighborhood at a certain field ridge are [[49.8, 50.1, 50.0], [49.9, 50.2, 50.1], [50.0, 50.3, 50.2]]. Therefore, Gx = (49.8 × (-1) + 50.1 × 0 + 50.0 × 1) + (49.9 × (-2) + 50.2 × 0 + 50.1). ×2)+(50.0×(-1)+50.3×0+50.2×1)=0.8, Gy=(49.8×(-1)+50.1×(-2)+50.0×(-1))+(49.9×0+50.2×0+50.1×0)+(50.0×1+50.3×2+50.2×1)=0.8, gradient magnitude G=√(0.8²+0.8²)=≈1.13, representing an elevation change of approximately 1.13 meters per meter at this location, indicating a region of drastic change, which is likely the boundary of the plot.

[0049] When generating a topographic gradient map, the gradient magnitude G is mapped to a grayscale value of 0-255 (formula: grayscale value = G × 255 / G_max, where G_max is the maximum gradient magnitude of the entire DEM, and in this example, G_max = 1.5). For example, the grayscale value corresponding to G = 1.13 is 1.13 × 255 / 1.5 ≈ 191. The higher the grayscale value, the greater the gradient. In the topographic gradient map, field ridges and plot edges appear as bright stripes, while flat farmland areas appear as dark tones, clearly distinguishing areas with drastic and gradual elevation changes.

[0050] The Canny edge detection algorithm is applied to the terrain gradient map to extract the contours of terrain abrupt changes and generate a preliminary edge contour map. Although the terrain gradient map highlights areas of elevation change, it still suffers from noise interference (such as individual high-gradient pixels) and edge discontinuities (such as breaks in the outline of field ridges). Therefore, the Canny edge detection algorithm is needed to further extract continuous and accurate terrain abrupt change contours. The Canny algorithm achieves high-quality edge extraction through multi-step processing. Core steps include Gaussian filtering, gradient calculation and direction determination, non-maximum suppression, double thresholding, and edge connection. Each step requires setting adaptation parameters based on the characteristics of the terrain gradient map.

[0051] First, perform Gaussian filtering. The purpose is to smooth the noise in the gradient map (such as isolated high-gradient pixels caused by DEM elevation errors). Select a 3×3 Gaussian kernel with σ = 1.0 (the Gaussian kernel function is G(x,y)=(1 / (2πσ²))×e^(-(x²+y²) / (2σ²))), and perform weighted averaging on each pixel of the gradient map. For example, if the gray value of a certain noise pixel is 200 (the gray values of surrounding pixels are all 50), after Gaussian filtering, the gray value of this pixel becomes (50×0.05 + 50×0.2 + 50×0.05 + 50×0.2 + 200×0.4 + 50×0.2 + 50×0.05 + 50×0.2 + 50×0.05) = 85. The noise is effectively suppressed, and at the same time, the gradient features of the main edges such as ridges are retained.

[0052] Secondly, recalculate the gradient direction based on the already smoothed gradient map (consistent with the calculation logic of the Sobel operator in the previous step), and discretize the gradient direction into four directions: 0°, 45°, 90°, and 135° (for subsequent non-maximum suppression). For example, if the gradient direction of a certain ridge pixel is 45°, it means that the elevation change direction of this pixel is southeast-northwest, which is consistent with the trend of the ridge.

[0053] Next, perform non-maximum suppression (NMS). The core is to refine the edges - for each pixel, judge whether it is a local maximum along the gradient direction. If so, retain it; otherwise, suppress it (set the gray value to 0) to avoid "thick lines" on the edges. For example, for a pixel with a gradient direction of 45° and a gray value of 191, the gray values of the two adjacent pixels along the 45° direction are 150 and 160 respectively, both of which are less than 191. Therefore, this pixel is retained; if the gray values of the adjacent pixels are 200 and 180 respectively, the current pixel is suppressed. Finally, the ridge edge is refined from a "bright band 3 - 5 pixels wide" to a "line 1 pixel wide".

[0054] Then, perform double-threshold processing. Set a high threshold Th and a low threshold Tl (determined according to the gray value distribution of the gradient map. In this example, Th = 0.7×G_max_gray = 0.7×255 = 178.5, Tl = 0.3×G_max_gray = 76.5, and G_max_gray is the maximum gray value 255 of the gradient map). Divide the pixels into three categories: pixels with gray value ≥ Th are "strong edges" (determined as boundaries), pixels with gray value < Tl are "non-edges" (discarded), and pixels with Tl ≤ gray value < Th are "weak edges" (need further judgment). For example, if the gray value of a certain ridge pixel is 190 ≥ 178.5, it is determined as a strong edge; if the gray value of a pixel in a flat area is 60 < 76.5, it is determined as a non-edge; if the gray value of a pixel in a transition area is 100, it is determined as a weak edge.

[0055] Finally, edge concatenation is performed, connecting weak edges with strong edges to form a continuous contour. If a weak edge is directly adjacent to a strong edge (within an 8-neighborhood), the weak edge is considered an edge and retained; if a weak edge is isolated and has no adjacent strong edges, it is discarded. For example, if a strong edge of a field ridge has a 2-pixel breakpoint, and a weak edge (grayscale value 100) at the breakpoint is adjacent to the strong edge, connecting them forms a continuous field ridge contour, avoiding edge breaks caused by low local gradients.

[0056] After the above processing, in the generated preliminary edge contour map, the boundaries of the plots (such as field ridges, the dividing line between plots and roads) appear as continuous white lines (grayscale value 255), and the background (flat farmland, internal area) is black (grayscale value 0). However, there may still be a small number of short and thin interference edges (such as the outline of small slopes), which need to be further optimized in subsequent steps.

[0057] The segmentation threshold is dynamically calculated based on local terrain features, and an adaptive threshold segmentation algorithm is used to transform the preliminary edge contour map into a binary boundary map. In the initial edge contour map, the edge gray values ​​of different regions differ (e.g., the edge gray values ​​are lower in flat areas and higher in sloping areas). If a fixed threshold is used for segmentation, problems such as "missing weak edges in flat areas" or "falsely detecting interfering edges in sloping areas" are likely to occur. The adaptive threshold segmentation algorithm solves the limitations of the fixed threshold by dynamically adjusting the threshold according to the gray-level characteristics of local areas. The core logic is to "use the gray-level statistics of the local area (such as the mean or median) as the segmentation threshold for that area," ensuring that the boundaries of different terrain regions can be accurately extracted.

[0058] First, determine the local window size. The window size needs to balance "representativeness of local features" and "computational efficiency". Based on the 1m×1m grid scale of the DEM, a 5×5 local window is selected (covering a 5m×5m area on the ground, which can reflect local terrain differences without making the threshold too smooth due to an excessively large window).

[0059] Then, the segmentation threshold for each local window is calculated using the "local mean method," where the segmentation threshold T(i,j) for each pixel is equal to the average gray value of all pixels in its 5×5 neighborhood. For example, in a flat area, the gray values ​​of edge pixels in the 5×5 neighborhood of a pixel are concentrated between 80 and 100, while the gray values ​​of background pixels are concentrated between 0 and 20. The neighborhood mean is approximately (80+90+100+0+20+...+0) / 25≈35, therefore the threshold T(i,j) for this pixel is 35. In a sloping area, the gray values ​​of edge pixels in the 5×5 neighborhood of a pixel are concentrated between 150 and 180, while the gray values ​​of background pixels are concentrated between 30 and 50. The neighborhood mean is approximately (150+160+180+30+50+...+40) / 25≈85, therefore the threshold T(i,j) for this pixel is 85.

[0060] Then, binarization is performed, comparing the grayscale value of each pixel with the corresponding local threshold T(i,j): if the pixel grayscale value > T(i,j), it is determined to be a "boundary pixel" and assigned a value of 1 (white, RGB value (255,255,255)); if the pixel grayscale value ≤ T(i,j), it is determined to be a "background pixel" and assigned a value of 0 (black, RGB value (0,0,0)). For example, in a flat area, a pixel with a grayscale value of 90 > 35 is assigned a value of 1 (preserving the boundary); in a sloping area, a pixel with a grayscale value of 70 ≤ 85 is assigned a value of 0 (suppressing interfering edges).

[0061] By using adaptive threshold segmentation, the binarized boundary map overcomes the limitations of fixed thresholds: weak edges (grayscale values ​​80-100) in flat areas are accurately preserved due to the low threshold (35), while interfering edges (grayscale values ​​60-70) in sloping areas are effectively removed due to the high threshold (85). In the final generated binarized boundary map, the main boundaries of the plots (such as field ridges and the boundary between plots and roads) appear as continuous white lines without obvious breaks or interfering lines, with a pure black background, laying the foundation for subsequent morphological processing.

[0062] Morphological closing operations are performed on the binary boundary map to fill the contour gaps and smooth the boundaries, generating an optimized binary boundary map. While binarized boundary maps remove most interference, two types of problems may still exist: first, "contour gaps" (such as 1-2 pixel wide gaps at the edge of a field ridge due to low local gradients); and second, "boundary burrs" (such as irregular protrusions of edge pixels). These problems can lead to discontinuous or irregular boundaries generated by subsequent vectorization, which need to be addressed through morphological closing operations. Morphological closing operations are a combination of "dilation followed by erosion." Dilation is used to fill gaps, and erosion is used to smooth boundaries. The combination of the two can both repair gaps and maintain the overall shape and position of the boundary.

[0063] First, an expansion operation is performed to fill the gaps at the boundaries and the edges of broken connections. A 3×3 rectangular structuring element is selected (the structuring element is a matrix of all 1s, i.e., [[1,1,1],[1,1,1],[1,1,1]]). The center of the structuring element is aligned with each pixel of the binarized image. If there is at least one "1" (boundary pixel) within the coverage area of ​​the structuring element, the center pixel is assigned a value of 1. For example, if there is a one-pixel-wide gap at the boundary of a field ridge (the gap pixel is 0, and the surrounding 4 pixels are 1), when the center of the structuring element is aligned with the gap pixel, because there are 4 "1"s within the coverage area, the gap pixel is assigned a value of 1, and the gap is filled; at the same time, the small depressions on the boundary are also filled, and the edge becomes more continuous.

[0064] The subsequent erosion operation aims to eliminate boundary burrs generated during the expansion process and restore the boundary to near its original width. A 3x3 rectangular structuring element, identical to the one used for expansion, is selected, with its center aligned with the pixel. The center pixel is assigned a value of 1 only if all pixels within the structuring element's coverage area are 1; otherwise, it is assigned a value of 0. For example, a 1-pixel-wide burr appearing on the boundary after expansion (the burr pixel is 1, and the surrounding 3 pixels are 0) is removed when the structuring element's center is aligned with the burr pixel because there are 0 pixels within its coverage area; the burr is thus removed. Simultaneously, the overall outline of the boundary does not shift significantly, and the direction and position of the field ridge remain unchanged.

[0065] The optimized boundary binary image after morphological closing operation has three significant characteristics: First, the boundary gaps are completely filled, and the main boundaries such as field ridges form complete closed contours without gaps of 1-2 pixels wide; second, boundary burrs are eliminated, edge smoothness is improved, and pixel-level irregular protrusions are reduced by more than 90%; third, the boundary width is uniform, basically maintaining a width of 1-2 pixels, which facilitates subsequent vectorization extraction. For example, before processing, a certain field ridge boundary had 3 gaps of 1 pixel and 5 burrs. After processing, all gaps were filled, and only 1 burr remained, significantly improving boundary continuity and smoothness.

[0066] The optimized binary boundary map is converted into an initial vector boundary map through vectorization and polygon fitting algorithms.

[0067] The optimized boundary binary map is in raster format (pixel matrix), which cannot be directly used for land area calculation and planning design. It needs to be converted into vector format (geometric objects composed of points, lines and surfaces). That is, the boundary lines are extracted by vectorization, and then the vector boundary of the land parcel is generated by polygon fitting.

[0068] First, vectorization is performed, the core of which is "edge tracing"—finding the first boundary pixel (a pixel with a grayscale value of 1) in the binary image as the starting point, recording its coordinates (e.g., (x0,y0)=(25,30), unit: pixels, which needs to be converted to geographic coordinates later). Then, following the 8-neighborhood order (top, top right, right, bottom right, bottom, bottom left, left, top left), the next adjacent boundary pixel is searched. If found, the coordinates are recorded and the current pixel is marked as "tracked". This process is repeated until the starting point is returned, forming a closed boundary line. For example, in the 8-neighborhood of the starting point (25,30), the right pixel (26,30) is the boundary pixel. The coordinates of this pixel are recorded and the pixel is moved to (26,30). The next boundary pixel in its 8-neighborhood is then searched. After continuous tracing, the coordinate sequence is obtained as: (25,30), (26,30), (27,31), (28,32),..., (25,30), forming a closed field ridge boundary line. For a binary map containing multiple plots, the edge tracing process needs to be repeated to extract the independent boundary lines of each plot and avoid boundary confusion.

[0069] Then, a coordinate transformation is performed to convert the pixel coordinates to geographic coordinates (consistent with the DEM's coordinate system, such as WGS84 latitude and longitude). The transformation formula is based on the DEM's grid scale and the geographic coordinates of the top-left corner: Assuming the geographic coordinates of the top-left corner pixel of the DEM are (Lon0, Lat0) = (118.5000°E, 30.2150°N), and the grid scale is 1m / pixel (i.e., each pixel corresponds to 1 meter of ground), then the geographic coordinates corresponding to a pixel coordinate (x, y) are: Longitude Lon = Lon0 + x × (1m / (111319.9m / °E)) (111319.9m is the approximate distance corresponding to 1 degree of longitude at the equator). Latitude Lat = Lat0 - y × (1m / (111319.9m / °N)) (Latitude is positive when it points north, and the pixel y-axis points downwards, so subtraction is used).

[0070] For example, the pixel coordinates (25, 30) correspond to the longitude = 118.5000°E + 25 × (1 / 111319.9)°E ≈ 118.5002°E, and the latitude = 30.2150°N - 30 × (1 / 111319.9)°N ≈ 30.2147. °N, completes the conversion from pixel coordinates to geographic coordinates.

[0071] Finally, polygon fitting is performed. The coordinate sequence obtained from edge tracking contains a large number of redundant points (such as continuous straight line segments containing dozens of pixels), which need to be simplified using the Douglas-Peucker algorithm to reduce the number of points while maintaining the boundary shape. The core of this algorithm is to set a "tolerance" (i.e., the maximum allowable deviation between the simplified boundary and the original boundary, set to 0.5m in accordance with land surveying accuracy requirements). The specific process is as follows: connect the starting point and the ending point of the coordinate sequence to form a straight line, calculate the distance from all intermediate points to this line, if the maximum distance > the tolerance, then retain the point with the largest distance as the dividing point, dividing the sequence into two segments; repeat this process until the distance from all intermediate points to the corresponding line is ≤ the tolerance, finally obtaining the simplified polygon vertices. For example, if the original coordinate sequence contains 50 points, after fitting, 8 key vertices are retained, and the deviation between the simplified boundary and the original boundary is ≤ 0.3m, meeting the accuracy requirements of land planning and design.

[0072] The simplified polygon vertices are connected sequentially to generate an initial vector boundary map. The vector format adopts the industry standard SHP format, which includes geometric information (polygon vertex coordinates) and attribute information (such as plot ID "LD-2025001", preliminary area "3333㎡", and boundary type "field ridge"), providing a vector data foundation for subsequent topology verification and area calculation.

[0073] S204, perform topological relationship verification and fragment patch fusion on the initial vector boundary map, use dynamic weight compensation algorithm to correct elevation error, and output the optimized land parcel vector boundary; Specifically, it can check the topological relationships of polygons in the initial vector boundary graph, verify closure and adjacency relationships, and generate a list of topological errors; The polygons in the initial vector boundary map are the core geometric representation of the land parcels. Their topological relationships (closure and adjacency) directly affect the accuracy of area calculation and subsequent planning and design. If topological errors exist (such as unclosed boundaries, overlapping adjacent parcels, or gaps), it will lead to deviations in area calculation or confusion in the spatial relationships of the parcels. Topological relationship verification needs to be carried out separately for "closure" and "adjacency," and verification thresholds should be set in conjunction with land survey accuracy requirements (boundary error ≤ 0.1 meters).

[0074] The core of closure verification is to ensure that the coordinates of the first and last vertices of the polygon coincide, avoiding boundary breaks due to missed vertex recording or coordinate deviation. The specific method is as follows: traverse the vertex coordinate sequence of each polygon (e.g., the vertex sequence of plot LD-2025001 is P1(118.5000°E,30.2150°N), P2(118.5200°E,30.2150°N), P3(118.5200°E,30.2000°N), P4(118.5000°E,30.2000°N), P5(118.5001°E,30.2149°N)), and calculate the Euclidean distance between the first and last vertices (P1 and P5). The distance calculation formula needs to be combined with the conversion relationship between latitude and longitude and meters (1° longitude corresponds to 111319.9 meters and 1° latitude corresponds to 111319.9 meters at the Earth's equator. Longitude distance needs to be multiplied by the cosine value of latitude for correction), that is: d=√[(Lon5-Lon1)²×(111319.9×cosLat_avg)²+(Lat5-Lat1)²×(111319.9)²], where Lat_avg is the average latitude of P1 and P5 (30.21495°), and cos30.21495°≈0.863. Substituting the coordinates, we get: Lon5-Lon1=0.0001°E, corresponding to a distance of 0.0001×111319.9×0.863≈9.61 meters; Lat5-Lat1=-0.0001°N, corresponding to a distance of 0.0001×111319.9≈11.13 meters; the total distance d≈√(9.61²+11.13²)≈14.7 meters, which far exceeds the preset closure threshold of 0.1 meters, and is judged as a "closure error".

[0075] Adjacency verification requires checking whether there are "overlaps" or "gaps" on the boundaries of adjacent polygons. Topological overlay analysis is used: the spatial intersection of the boundary segments of two adjacent plots is calculated. If the intersection length is more than 10% of the length of a certain boundary segment (e.g., the intersection length of boundary segment AB of plot A and boundary segment CD of plot B is 2 meters, and the total length of AB is 15 meters, accounting for 13.3%), it is judged as "overlap error"; if the minimum distance between the two boundary segments exceeds 0.2 meters (e.g., the distance between the boundary segments of plot C and plot D is 0.3 meters, and there is a blank area), it is judged as "gaps error".

[0076] The final generated topology error list must include "Error ID, Error Type, Affected Plot ID, Error Location Coordinates, and Error Quantification Value", for example: "Error ID: TOP-001, Type: Closure Error, Plot ID: LD-2025001, Location: P1(118.5000°E,30.2150°N)-P5(118.5001°E,30.2149°N), End-to-End Distance: 14.7 meters; Error ID: TOP-002, Type: Overlap Error, Plot ID: LD-2025001 / LD-2025002, Location: AB segment (118.5050°E,30.2100°N)-(118.5100°E,30.2100°N), Overlap Length: 2 meters", providing a clear basis for subsequent corrections.

[0077] Based on the list of topology errors, perform boundary trimming and node adjustment to eliminate overlaps and gaps, and generate a topologically correct vector boundary map; Topology error correction should follow the "principle of minimum adjustment" to avoid excessive modifications that could cause boundaries to deviate from the actual terrain. Targeted methods should be used for different error types. For closure errors, if the distance between the first and last points is small (0.1-0.5 meters), the coordinates of the last point are directly corrected to the coordinates of the first point. For example, if the distance between the first and last points of a polygon is 0.3 meters, the last point P5 (118.5001°E, 30.2149°N) is adjusted to the first point P1 (118.5000°E, 30.2150°N), reducing the closure error to 0.02 meters. If the distance is large (e.g., 14.7 meters), it is necessary to check whether there are any missing or out-of-order points in the vertex sequence. If the missing vertex P6 (118.5000°E, 30.2000°N) is found, the vertex sequence after supplementation is P1-P2-P3-P4-P6-P1, reducing the distance between the first and last points to 0.01 meters, which meets the closure requirements.

[0078] For overlap errors, "boundary trimming" is used: extract the coordinates of the overlapping boundary segments (e.g., AB segment of plot A: A(118.5050°E, 30.2100°N), B(118.5100°E, 30.2100°N); CD segment of plot B: C(118.5060°E, 30.2100°N), D(118.5110°E, 30.2100°N), and then...). (Segment CB), calculate the midpoint coordinates of the overlapping segment ((118.5060+118.5100) / 2°E,30.2100°N)=(118.5080°E,30.2100°N), correct point B of plot A to the midpoint, correct point C of plot B to the midpoint, after trimming segment AB becomes A-midpoint, segment CD becomes midpoint-D, the overlap is completely eliminated, and the boundary position conforms to the original terrain direction.

[0079] For gap errors, "node snapping" is used: find the nearest boundary nodes on both sides of the gap (such as node E (118.5150°E, 30.2050°N) for plot C and node F (118.5152°E, 30.2051°N) for plot D, with a spacing of 0.3 meters), calculate the midpoint of the line connecting the two points (118.5151°E, 30.20505°N), correct E and F to the midpoint coordinates, close the gap, and the continuity of adjacent boundaries is not affected.

[0080] After correction, the topological relationships need to be re-verified to ensure that all errors are eliminated—for example, the closure error of plot LD-2025001 is ≤0.02 meters, adjacent plots have no overlap or gaps, and a topologically correct vector boundary map is generated to lay the foundation for subsequent fragment patch fusion.

[0081] Calculate the area and perimeter of each patch in the vector boundary map, identify and fuse fragmented patches based on the area threshold, and generate a fused vector boundary map; The area and perimeter of a plot are the core criteria for determining whether it is a "fragmented plot". Fragmented plots (redundant plots with too small an area) will interfere with land planning and design and need to be integrated into the adjacent main plots.

[0082] The area calculation uses the shoelace formula, which is applicable to solving the area of ​​any polygon. The formula is: S=0.5×|Σ(Lon_i×Lat_{i+1}-Lon_{i+1}×Lat_i)| (i=1 to n, Lon_{n+1}=Lon_1, Lat_{n+1}=Lat_1). The result needs to be converted to square meters using the latitude-longitude conversion factor (conversion factor = 111319.9²×cosLat_avg, where cosLat_avg is the average latitude cosine of the polygon vertices). For example, given the vertex coordinates of a certain polygon as Q1(118.5000°E, 30.2150°N), Q2(118.5010°E, 30.2150°N), Q3(118.5010°E, 30.2140°N), and Q4(118.5000°E, 30.2140°N), substituting these coordinates into the formula, we can calculate: Σ(Lon_i×Lat_{i+1})=118.5000×30.2150+118.5010×30.2140+118.5010×30.2140+1 18.5000×30.2150≈14371.4;Σ(Lon_{i+1}×Lat_i)=30.2150×118.5010+30.2150×118.5010+30.2140×118.5000+30.2140×118.5000≈14371.2;S=0.5×|14371.4-14371.2|×(111319.9²×0.863)≈980 square meters (error ≤2%, close to the theoretical value of 1000 square meters).

[0083] The perimeter is calculated as the sum of the Euclidean distances between adjacent vertices. For example, the distance between Q1 and Q2 is approximately 0.001°E × 111319.9 × 0.863 ≈ 96.1 meters, Q2-Q3 is approximately 0.001°N × 111319.9 ≈ 111.3 meters, Q3-Q4 is approximately 96.1 meters, Q4-Q1 is approximately 111.3 meters, and the total perimeter is approximately 414.8 meters.

[0084] Fragmented land parcels are identified based on an area threshold, which needs to be set in conjunction with the minimum land use unit in land planning. The minimum plot area threshold in rural land planning is typically 66.7 square meters (0.1 mu), while for urban construction land it is 100 square meters. In this example, 66.7 square meters is used. All land parcels are iterated over. If a land parcel has an area of ​​60 square meters (less than 66.7 square meters), it is identified as a fragmented land parcel, and its adjacent land parcel is plot LD-2025002 with an area of ​​3000 square meters.

[0085] The fusion method is "boundary merging": extract the common boundary between fragmented patches and adjacent patches (if there is no boundary, take the shortest distance boundary), delete the vertices of the common boundary, and merge the vertex sequence of the fragmented patch into the vertex sequence of the adjacent patch. For example, after merging the vertices Q1-Q2-Q3-Q4 of fragmented patch with the vertices R1-R2-Q2-Q3-R3-R4 of plot LD-2025002, the new vertex sequence is R1-R2-Q1-Q4-Q3-R3-R4. The area of ​​the merged patch is 3000 + 60 = 3060 square meters, the boundary is continuous without breaks, and a fused vector boundary map is generated, eliminating the interference of redundant small plots on subsequent processing.

[0086] Using a dynamic weight compensation algorithm, the vector boundary points are adjusted and corrected based on the elevation values ​​of the digital elevation model to generate the elevation-corrected vector boundary. The core logic of the dynamic weight compensation algorithm is as follows: In the digital elevation model (DEM), regions with small elevation variance (such as flat farmland) have high elevation accuracy at their corresponding boundary points and are assigned higher weights; regions with large variance (such as steep slopes and gullies) have low elevation accuracy and are assigned lower weights. By weighted adjustment, the influence of high-error regions on the boundary position is reduced, thereby improving the spatial accuracy of the boundary.

[0087] Specific implementation steps: Obtain DEM elevation and variance: Each vector boundary point corresponds to a grid in the DEM. Extract the elevation values ​​of this grid and its 3×3 neighboring grids, and calculate the elevation variance (the smaller the variance, the higher the elevation accuracy). For example, boundary point P (118.5050°E, 30.2100°N) falls within DEM grid G. The elevation values ​​of grid G ​​and its neighboring grids are [50.1, 50.2, 50.3, 50.2, 50.2, 50.1, 50.3, 50.2, 50.2], with an average elevation Havg = 50.2 meters. The variance σ² = Σ(Hi - Havg)² / 9 ≈ (0.1²×4 + 0.1²×2) / 9 ≈ 0.0067.

[0088] Calculate the weights: The weights are inversely proportional to the variance, and the formula is w=1 / σ², where w=1 / 0.0067≈149.3 (the smaller the variance, the larger the weight, and the higher the priority of boundary point adjustment).

[0089] Weighted adjustment correction: For each boundary point and its two adjacent boundary points (preceding point P_prev, subsequent point P_next), the plane coordinates (Lon, Lat) are corrected. The adjustment formula is as follows: Lon_corr=(w×Lon+w_prev×Lon_prev+w_next×Lon_next) / (w+w_prev+w_next); Lat_corr=(w×Lat+w_prev×Lat_prev+w_next×Lat_next) / (w+w_prev+w_next); Where w_prev and w_next are the weights of P_prev and P_next, respectively (e.g., the variance of P_prev is σ_prev²=0.01, w_prev=100; the variance of P_next is σ_next²=0.005, w_next=200). Substituting the original coordinates of P (Lon=118.5050°E, Lat=30.2100°N), P_prev (Lon=118.5040°E), and P_next (Lon=118.5060°E), we calculate: Lon_corr=(149.3×118.5050+100×118.5040+200×118.5060) / (149.3+100+ 200)≈118.5052°E; Lat_corr≈30.2101°N.

[0090] After correction, the elevation error of boundary point P decreased from ±0.15 meters to ±0.08 meters, which is closer to the actual terrain. After all boundary points are corrected, a vector boundary with corrected elevation is generated, providing a high-precision coordinate basis for subsequent smoothing processing.

[0091] The vector boundaries of the elevation correction are smoothed using Bézier curves, and the optimized parcel vector boundaries are output.

[0092] The vector boundary after elevation correction is still a "broken line" shape with obvious sharp corners (such as the corners of field ridges and the turning points of plot edges), which does not conform to the continuous shape of natural terrain. It needs to be smoothed by Bézier curves to eliminate sharp corners and maintain the spatial accuracy of the boundary.

[0093] A cubic Bézier curve is used, which is defined by four control points (start point P0, end point P3, and intermediate control points P1 and P2). The curve equation is: B(t)=(1-t)³P0+3(1-t)²tP1+3(1-t)t²P2+t³P3 (t∈[0,1]). The smoothness is controlled by adjusting the position of the intermediate control points to ensure that the curve follows the original boundary direction and the area error is ≤0.5%.

[0094] Specific implementation: Grouping and Control Point Selection: The elevation-corrected boundary vertex sequence is divided into groups of three consecutive vertices (e.g., P0-P1-P2, P2-P3-P4, etc.), and each group generates a Bézier curve. Intermediate control points follow the principle of "consistent tangent direction": For vertex P1 (connecting P0-P1-P2), the preceding control point P1_prev is the extension of the line connecting P1 and P0, with a distance from P1 equal to 1 / 3 of the length of P0-P1; the following control point P1_next is the extension of the line connecting P1 and P2, with a distance from P1 equal to 1 / 3 of the length of P1-P2. For example, P0 (118.5000°E, 30.2150°N), P1 (118.5050°E, 30.2150°N), P2 (118.5050°E, 30.2100°N), the length of P0-P1 is approximately 96.1 meters, P1_prev=P1+(P0-P1) / 3≈(118.5033°E, 30.2150°N); the length of P1-P2 is approximately 111.3 meters, P1_next=P1+(P2-P1) / 3≈(118.5050°E, 30.2133°N).

[0095] Curve generation and verification: Substituting into the curve equation, a smooth curve from P0 to P2 is generated when t changes from 0 to 1. For example, when t=0.5, B(0.5)≈(118.5039°E,30.2131°N), and this point is located inside the original broken line, forming a natural corner. Verification of area error after smoothing: The original plot area was 3060 square meters, and after smoothing, it became 3058 square meters, with an error of 0.06%, meeting the accuracy requirements.

[0096] The final output is the optimized parcel vector boundary, which is a continuous and smooth curve with no topological errors, fragmented patches, or elevation deviations. It can be directly used for subsequent coordinate registration with the geographic information system base map.

[0097] S205, the optimized land parcel vector boundary is registered with the geographic information system base map, the projected area is calculated through integration, and a final measurement report including the terrain undulation coefficient is generated.

[0098] Specifically, the optimized land parcel vector boundary can be transformed to the coordinate system of the geographic information system base map, and a seven-parameter coordinate transformation model can be used for projection transformation to generate the vector boundary after coordinate transformation. The optimized plot vector boundary defaults to the native coordinate system of the UAV-collected data—the WGS84 geodetic coordinate system (latitude and longitude coordinates, unit: degrees). However, the Geographic Information System (GIS) base map typically uses a Cartesian coordinate system suitable for local planning (such as the Beijing 54 coordinate system or the Xi'an 80 coordinate system). Therefore, a seven-parameter coordinate transformation model is needed to achieve accurate cross-coordinate system transformation. The seven-parameter model is the core model for geodetic coordinate system transformation, containing seven key parameters: three translation parameters (ΔX, ΔY, ΔZ, representing the offset of the origin of the two coordinate systems in the spatial Cartesian coordinate system), three rotation parameters (εx, εy, εz, representing the small rotation angles of the two coordinate systems around the X, Y, and Z axes, unit: seconds), and one scale parameter (m, representing the scale difference coefficient between the two coordinate systems). All parameters need to be obtained through calibration using known control points (features that simultaneously possess coordinates in both the WGS84 and base map coordinate systems).

[0099] Taking the conversion from "WGS84 coordinate system to Xi'an 80 coordinate system (3-degree zone, central meridian 117°E)" as an example, the seven parameters are first calculated using three campus control points (such as the southeast corner of the playground and the northwest corner of the teaching building): Translation parameters: ΔX = 6.2 meters (X-axis offset, i.e., the distance between the origins of the two coordinate systems on the X-axis), ΔY = -4.5 meters (Y-axis offset), ΔZ = 2.1 meters (Z-axis offset); Rotation parameters: εx = 1.8 seconds (rotation around the X-axis, positive values ​​are counterclockwise), εy = -2.3 seconds (rotation around the Y-axis), εz = 1.1 seconds (rotation around the Z-axis); Scale parameter: m = 1 + 1.5 × 10 -6 (i.e., 1.0000015, representing the scale magnification factor of the Xi'an 80 coordinate system relative to WGS84, the value is close to 1, indicating that the scale difference is extremely small).

[0100] The conversion process consists of two steps: Geodetic coordinates to spatial rectangular coordinates: The optimized boundary points' WGS84 latitude and longitude (e.g., point P: Lon=118.5000°E, Lat=30.2150°N, H=50.2 meters, where H is the geodetic height) are converted to WGS84 spatial rectangular coordinates (X, Y, Z, unit: meters). The formula is based on the Earth ellipsoid parameters (WGS84 ellipsoid semi-major axis a=6378137 meters, flattening f=1 / 298.257223563), and the calculated values ​​are X=3345678.12 meters, Y=118567.34 meters, and Z=2890123.56 meters. Seven-parameter transformation: Substituting the seven-parameter transformation formula, the WGS84 spatial rectangular coordinates are converted to Xi'an 80 spatial rectangular coordinates. The formula is as follows: X'=X+ΔX+εy×Z-εz×Y+m×X; Y'=Y+ΔY+εz×X-εx×Z+m×Y; Z'=Z+ΔZ+εx×Y-εy×X+m×Z.

[0101] The rotation parameters need to be converted to radians first (1 second ≈ 4.848 × 10⁻⁶). -6 (radians), such as εx = 1.8 seconds ≈ 8.726 × 10⁻⁶ seconds. -6 In radians. Substituting the coordinates of point P, we get X' = 3345678.12 + 6.2 + (-2.3 × 4.848e-6) × 2890123.56 - (1.1 × 4.848e-6) × 118567.34 + 1.0000015 × 3345678.12 ≈3345684.35 meters, Y'≈118562.86 meters; Spatial rectangular coordinates to plane rectangular coordinates: The spatial rectangular coordinates of Xi'an 80 are converted to plane rectangular coordinates through Gauss projection (3-degree zone, zone number 39, central meridian 117°E). The final plane coordinates of point P in Xi'an 80 are (X=3345684.35 meters, Y=39118562.86 meters, with zone number 39 added before the Y coordinate).

[0102] The transformation is completed by traversing all boundary points, generating a vector boundary after coordinate transformation. The boundary coordinate accuracy is ≤0.1 meters, which meets the requirements for GIS base map overlay.

[0103] Select control points in the base map and perform least-squares registration with the vector boundary after coordinate transformation to generate the registered vector boundary. The vector boundary after coordinate transformation may have slight offsets (usually ≤0.5 meters) due to seven-parameter errors and UAV positioning deviations. The offsets need to be further eliminated by least squares registration. The core is to use high-precision control points (such as road intersections and fixed building corners with coordinate accuracy ≤0.05 meters) on the GIS base map to establish error equations, solve for the optimal transformation parameters, and align the boundary with the features on the base map.

[0104] 1. Control point selection and data acquisition Select at least 3 evenly distributed control points from the GIS base map (the number must be greater than or equal to the number of transformation parameters; planar registration requires 4 parameters, so 4 are selected). Record the "base map coordinates (X-base, Y-base)" and "coordinates of the corresponding point on the transformed boundary (X-rotation, Y-rotation)" for each control point. Example data: Control point C1: Base map coordinates (3345000.00 meters, 39118000.00 meters), converted coordinates (3345002.10 meters, 39117998.50 meters). Control point C2: Base map coordinates (3345500.00 meters, 39118500.00 meters), transformed coordinates (3345501.80 meters, 39118502.30 meters). Control point C3: Base map coordinates (3346000.00 meters, 39118000.00 meters), converted coordinates (3346003.20 meters, 39117999.10 meters). Control point C4: Base map coordinates (3345500.00 meters, 39117500.00 meters), converted coordinates (3345502.50 meters, 39117498.80 meters).

[0105] 2. Establishment of the least squares registration model Planar registration uses a "similar transformation model", which includes four parameters: X-direction translation a, Y-direction translation b, rotation angle θ, and scale factor k. The transformation formulas are: X base = a + k × (X rotation × cosθ - Y rotation × sinθ); Y base = b + k × (X rotation × sinθ + Y rotation × cosθ); Because the offset is small, θ≈0 (cosθ≈1, sinθ≈θ, θ unit: radians), k≈1, the model can be simplified to a linear form, establishing the error equations: vxi=a+X_rotation×θ-Y_rotation×k'-X_base (k'=k-1, representing the small scale deviation); vyi=b+Y_rotation×θ+X_rotation×k'-Y_base. Here, vxi and vyi are the residuals in the X and Y directions (i.e., registration errors), and the goal is to minimize the sum of squared residuals Σ(vxi²+vyi²).

[0106] 3. Parameter Solving and Boundary Correction Parameters a, b, θ, and k' are solved using matrix operations. Substituting these parameters into the example data, we obtain: a = -2.08 meters (X-direction correction), b = 1.95 meters (Y-direction correction), θ ≈ -0.00012 radians (approximately -0.0069°, slight rotation), k' ≈ -1.2 × 10⁻⁶. -6 (Scale Correction). These parameters are used to correct all transformed boundary points. For example, the corrected coordinates of point P (transformed coordinates 3345684.35 meters, 39118562.86 meters) are: X registration = -2.08 + 1.0000012 × (3345684.35 × 1 - 39118562.86 × (-0.00012)) ≈ 3,345,682.26 meters; Y-registration = 1.95 + 1.0000012 × (3345684.35 × (-0.00012) + 39118562.86 × 1) ≈ 3,345,684.82 meters; After correction, the residuals of all control points are ≤0.08 meters, the registration accuracy meets the land planning requirements, and the registered vector boundary is generated.

[0107] Perform Gaussian projection area integration on the registered vector boundary to calculate the projected area of ​​the plot and generate the area calculation result; The registered vector boundary is located in a Cartesian coordinate system (Xi'an 80). The area of ​​the plot needs to be calculated by integrating the area of ​​the Gaussian projection. The Gaussian projection is an equal-angle transverse cylindrical projection. After the ellipsoidal plot is projected onto the plane, the area needs to be corrected for projection distortion by integration to ensure that the deviation from the actual ellipsoidal area is ≤0.1%.

[0108] 1. Gaussian projection zoning and scale factor The target site is located in the 3-degree zone of the Xi'an 80 coordinate system (zone number 39, central meridian 117°E). The scale factor m of the Gaussian projection varies with longitude (m=1 at the central meridian, and m increases with distance from the central meridian). Therefore, the scale factor at the center of the site needs to be calculated first. The longitude of the site center is Lon=118.5000°E, and the difference from the central meridian is ΔLon=1.5°. The scale factor formula is m=1+(ΔLon×π / 180)²×(a²×cos²Lat) / (2×N²), where Lat=30.2150°N (latitude of the site center), a=6378140 meters (semi-major axis of the Xi'an 80 ellipsoid), and N=a / √(1-f×sin²Lat) (radius of curvature of the trochanteric and vertex axes). The calculated value is m≈1.00003 (minimal projection distortion).

[0109] 2. Calculation of the area of ​​polygon integrals The registered vector boundary is treated as a closed polygon, and the area is calculated using the "trapezoidal integral method": the vertices of the polygon are arranged in clockwise order (e.g., A(XA,YA), B(XB,YB), C(XC,YC), D(XD,YD),A), and the polygon is divided into several trapezoids along the X-axis. The area of ​​each trapezoid is (upper base + lower base) × height / 2, and the total area is the sum of the areas of all trapezoids.

[0110] Example vertex coordinates (simplified): A (3345000.00 m, 39118000.00 m), B (3345500.00 m, 39118500.00 m), C (3346000.00 m, 39118000.00 m), D (3345500.00 m, 39117500.00 m) Integral calculation process: Divide the area along the X-axis from the minimum value XA = 3345000.00 meters to the maximum value XC = 3346000.00 meters, according to the X-coordinate of the vertex: Segment 1 (XA→XB): X changes from 3345000 to 3345500 meters. The left boundary Y = YA + (YB-YA) × (X-XA) / (XB-XA), and the right boundary Y = YD + (YC-YD) × (X-XA) / (XC-XA). Integrating, we get the area S1 = (YA+YB) × (XB-XA) / 2 = (118000+118500) × 500 / 2 = 59125000 square meters (Note: Y coordinate omits zone 39 for simplified calculation). Segment 2 (XB→XC): X changes from 3345500 to 3346000 meters. The left boundary Y = YB + (YC-YB) × (X-XB) / (XC-XB), and the right boundary Y = YD + (YC-YD) × (X-XB) / (XC-XB). Integrating, we get the area S2 = (YB+YC) × (XC-XB) / 2 = (118500 + 118000) × 500 / 2 = 59125000 square meters. The total area S_plane = S1 + S2 = 118,250,000 square meters (actual calculations require coordinate simplification and adjustment; the actual area is approximately 250,000 square meters, or 75 mu).

[0111] 3. Area Correction and Result Generation Due to scale distortion in Gaussian projection, the plane integral area needs to be corrected to the ellipsoidal area, using the formula S_ellipsoid = S_plane / m. Substituting m = 1.00003, we get S_ellipsoid = 250000 / 1.00003 ≈ 249992.5 square meters. The generated area calculation result is: "Gaussian projection area of ​​the plot: 250000.00 square meters (equivalent to 75.00 mu), actual ellipsoidal area: 249992.50 square meters, area error: 0.003%, meeting the accuracy requirements."

[0112] The terrain relief coefficient is calculated based on the digital elevation model, and the area calculation results and terrain parameters including at least the terrain relief coefficient are integrated to generate a final measurement report containing the terrain parameters.

[0113] The topographic relief coefficient is a core indicator that characterizes the flatness of a plot of land. It needs to be calculated in conjunction with a digital elevation model (DEM), and auxiliary topographic parameters such as average elevation and average slope should be extracted and integrated with the area results into a final report to provide complete data support for land planning.

[0114] 1. Calculation of terrain relief coefficient The topographic relief coefficient K = σ / H_avg, where σ is the standard deviation of the elevation of the plot area in the DEM (reflecting the degree of elevation dispersion), and H_avg is the average elevation (reflecting the overall elevation of the plot).

[0115] Extract DEM data: Select 100 evenly distributed DEM grid points (1m × 1m resolution) within the plot area, with elevation values ​​ranging from 48.5 to 52.3 meters. Example data: 48.5, 49.0, 49.2, 49.5, 50.0, 50.2, 50.5, 51.0, 51.5, 52.3 (10 points in total, for simplified calculation). Calculate the average elevation H_avg = (48.5 + 49.0 + 49.2 + 49.5 + 50.0 + 50.2 + 50.5 + 51.0 + 51.5 + ... 52.3) / 10≈50.17 meters; Calculate the standard deviation σ: First, calculate the sum of squares of the deviations of each elevation from the mean Σ(H_i-H_avg)²=(48.5-50.17)²+(49.0-50.17)²+...+(52.3-50.17)²≈12.79, variance=12.79 / (10-1)≈1.42, σ=√1.42≈1.19 meters; The topographic relief coefficient K = 1.19 / 50.17 ≈ 0.024 (K < 0.05 indicates slight undulation, 0.05-0.1 indicates moderate undulation, and > 0.1 indicates strong undulation. In this example, it is slight undulation, which is suitable for agricultural planting or construction land).

[0116] 2. Calculation of auxiliary terrain parameters Average slope: The slope of each grid point is calculated using the DEM (the elevation gradient is calculated using the Sobel operator, and the slope angle = arctan(gradient magnitude)). The slope range of 100 points is 0.5°-3.2°, and the average slope is ≈1.8°. Slope aspect distribution: According to the slope direction statistics, north-facing (0°-90°) accounts for 32%, south-facing (180°-270°) accounts for 28%, east-facing (90°-180°) accounts for 20%, and west-facing (270°-360°) accounts for 20%. South-facing slopes receive ample sunlight and are suitable for crop growth.

[0117] 3. Final measurement report generation The report uses a standardized format and includes six core modules: Project Overview: Name of the surveyed plot (East Plot of XX Village), Location (118.5000°-118.5200°E, 30.2000°-30.2150°N), Purpose of Survey (Agricultural Planting Planning), Data Source (DJI M300 drone + MicaSenseAltum sensor, DEM resolution 1m×1m). Area calculation results: Gaussian projection area 250,000.00 square meters (75.00 mu), ellipsoidal surface area 249,992.50 square meters, area error 0.003%, accuracy level 1; Topographic parameters: average elevation 50.17 meters, topographic relief coefficient 0.024 (slight undulation), average slope 1.8°, with north-south orientation predominant (60%). Conclusions and Recommendations: The plot has a flat terrain and good drainage, making it suitable for planting grain crops such as rice and wheat; the area accuracy meets the requirements of land planning and design and can be used as the basis for land use approval. Data attachments: Registered vector boundary map (SHP format), DEM grayscale image, control point coordinate table; Measurement unit and date: XX Surveying and Mapping Co., Ltd., October 1, 2025.

[0118] The report is also uploaded to the GIS management platform for use by planners.

[0119] Another embodiment of the present invention provides a land area measurement system for land planning and design, see [link to relevant documentation]. Figure 3 The system may include: The acquisition module 301 is used to acquire a sequence of aerial images of the target site using a drone equipped with a multispectral sensor, and to simultaneously record the positioning and attitude data of each frame of the image. Processing module 302 is used to perform edge sharpening processing and multi-view three-dimensional reconstruction on the sequence of aerial images to generate a digital elevation model with texture features. Analysis module 303 is used to perform terrain gradient analysis and edge contour extraction based on the digital elevation model, identify land parcel boundaries and generate an initial vector boundary map through an adaptive threshold segmentation algorithm; The verification module 304 is used to perform topological relationship verification and fragment patch fusion on the initial vector boundary map, correct elevation errors using a dynamic weight compensation algorithm, and output the optimized land parcel vector boundary. The registration module 305 is used to register the optimized land parcel vector boundary with the geographic information system base map, calculate the projected area through integration, and generate a final measurement report containing the terrain undulation coefficient.

[0120] The above description, based on the embodiments shown in the figures, details the structure, features, and effects of the present invention. The above description is only a preferred embodiment of the present invention, but the present invention is not limited to the scope of implementation shown in the figures. Any changes made in accordance with the concept of the present invention, or equivalent embodiments modified to have equivalent changes, that do not exceed the spirit covered by the specification and figures, should be within the protection scope of the present invention.

Claims

1. A method for measuring land area for land planning and design, characterized in that, The method includes: A sequence of aerial images of the target site is collected by a drone equipped with a multispectral sensor, and the positioning and attitude data of each frame of the image are recorded simultaneously. The sequence of aerial images is subjected to edge sharpening and multi-view 3D reconstruction to generate a digital elevation model with texture features. Based on the digital elevation model, terrain gradient analysis and edge contour extraction are performed, and the plot boundaries are identified and an initial vector boundary map is generated by an adaptive threshold segmentation algorithm. The initial vector boundary map is subjected to topological relationship verification and fragment patch fusion. The elevation error is corrected by a dynamic weight compensation algorithm, and the optimized land parcel vector boundary is output. The optimized land parcel vector boundary is registered with the geographic information system base map, the projected area is calculated through integration, and a final measurement report including the terrain undulation coefficient is generated.

2. The method according to claim 1, characterized in that, The process of acquiring a sequence of aerial images of the target site using a drone equipped with a multispectral sensor and simultaneously recording the positioning and attitude data of each frame includes: Plan the UAV flight path based on the scope and shape of the target site, and determine the flight path plan and sensor parameter configuration; Control the drone to fly according to the flight path plan, collect a sequence of aerial images through multispectral sensors, and record high-precision positioning and attitude data of each frame of the image to generate the original aerial images and positioning and attitude dataset; The original aerial images and positioning and attitude determination datasets are time-synchronized and quality-checked, blurry or missing data frames are removed, and the verified sequence of aerial images and synchronized positioning and attitude determination data are generated.

3. The method according to claim 2, characterized in that, The step of performing edge sharpening processing and multi-view 3D reconstruction on the sequence of aerial images to generate a digital elevation model with texture features includes: Radiometric correction and noise filtering preprocessing are performed on the verified sequence of aerial images to eliminate illumination differences and sensor noise, generating a preprocessed image sequence. The Laplacian edge sharpening algorithm is applied to the preprocessed image sequence to perform convolution operation, thereby enhancing the boundary features of ground objects and generating an edge-sharpened image sequence. Based on synchronous positioning and attitude determination data, a scale-invariant feature transformation algorithm is used to match key points in multi-view images and generate a set of matching feature points. A multi-view stereo vision algorithm is used to reconstruct a dense 3D point cloud from the matching feature point set, generating a preliminary 3D point cloud model. The texture information of the edge-sharpened image sequence is mapped onto the preliminary 3D point cloud model, and a digital elevation model with texture features is generated through surface reconstruction algorithms.

4. The method according to claim 3, characterized in that, The process of performing terrain gradient analysis and edge contour extraction based on the digital elevation model, identifying land parcel boundaries and generating an initial vector boundary map using an adaptive threshold segmentation algorithm includes: The terrain gradient of the digital elevation model is calculated, and the spatial derivative of the elevation value is solved using the Sobel operator to generate a terrain gradient map. The Canny edge detection algorithm is applied to the terrain gradient map to extract the contours of terrain abrupt changes and generate a preliminary edge contour map. The segmentation threshold is dynamically calculated based on local terrain features, and an adaptive threshold segmentation algorithm is used to transform the preliminary edge contour map into a binary boundary map. Morphological closing operations are performed on the binary boundary map to fill the contour gaps and smooth the boundaries, generating an optimized binary boundary map. The optimized binary boundary map is converted into an initial vector boundary map through vectorization and polygon fitting algorithms.

5. The method according to claim 4, characterized in that, The process of performing topological relationship verification and fragment patch fusion on the initial vector boundary map, correcting elevation errors using a dynamic weight compensation algorithm, and outputting optimized land parcel vector boundaries includes: Check the topological relationships of polygons in the initial vector boundary graph, verify closure and adjacency, and generate a list of topological errors; Based on the list of topology errors, perform boundary trimming and node adjustment to eliminate overlaps and gaps, and generate a topologically correct vector boundary map; Calculate the area and perimeter of each patch in the vector boundary map, identify and fuse fragmented patches based on the area threshold, and generate a fused vector boundary map; Using a dynamic weight compensation algorithm, the vector boundary points are adjusted and corrected based on the elevation values ​​of the digital elevation model to generate the elevation-corrected vector boundary. The vector boundaries of the elevation correction are smoothed using Bézier curves, and the optimized parcel vector boundaries are output.

6. The method according to claim 5, characterized in that, The step of registering the optimized land parcel vector boundary with the geographic information system base map, calculating the projected area through integration, and generating a final measurement report including terrain relief coefficients includes: The optimized land parcel vector boundary is transformed to the coordinate system of the geographic information system base map. A seven-parameter coordinate transformation model is used for projection transformation to generate the vector boundary after coordinate transformation. Select control points in the base map and perform least-squares registration with the vector boundary after coordinate transformation to generate the registered vector boundary. Perform Gaussian projection area integration on the registered vector boundary to calculate the projected area of ​​the plot and generate the area calculation result; The terrain relief coefficient is calculated based on the digital elevation model, and the area calculation results and terrain parameters including at least the terrain relief coefficient are integrated to generate a final measurement report containing the terrain parameters.

7. A land area measurement system for land planning and design, characterized in that, The system includes: The acquisition module is used to acquire a sequence of aerial images of the target site using a drone equipped with a multispectral sensor, and to simultaneously record the positioning and attitude data of each frame of the image. The processing module is used to perform edge sharpening and multi-view 3D reconstruction on the sequence of aerial images to generate a digital elevation model with texture features. The analysis module is used to perform terrain gradient analysis and edge contour extraction based on the digital elevation model, identify land parcel boundaries and generate an initial vector boundary map through an adaptive threshold segmentation algorithm; The verification module is used to perform topological relationship verification and fragment patch fusion on the initial vector boundary map, correct elevation errors using a dynamic weight compensation algorithm, and output the optimized land parcel vector boundary. The registration module is used to register the optimized land parcel vector boundary with the geographic information system base map, calculate the projected area through integration, and generate a final measurement report including the terrain undulation coefficient.

8. The system according to claim 7, characterized in that, The process of acquiring a sequence of aerial images of the target site using a drone equipped with a multispectral sensor and simultaneously recording the positioning and attitude data of each frame includes: Plan the UAV flight path based on the scope and shape of the target site, and determine the flight path plan and sensor parameter configuration; Control the drone to fly according to the flight path plan, collect a sequence of aerial images through multispectral sensors, and record high-precision positioning and attitude data of each frame of the image to generate the original aerial images and positioning and attitude dataset; The original aerial images and positioning and attitude determination datasets are time-synchronized and quality-checked, blurry or missing data frames are removed, and the verified sequence of aerial images and synchronized positioning and attitude determination data are generated.

9. A storage medium, characterized in that, The storage medium stores a computer program, wherein the computer program is configured to execute the method of any one of claims 1-6 when it is run.

10. An electronic device comprising a memory and a processor, characterized in that, The memory stores a computer program, and the processor is configured to run the computer program to perform the method of any one of claims 1-6.

Images