A Land Boundary Measurement Method Based on UAV Imagery

Abstract

This invention relates to the field of photogrammetry, specifically to a land boundary measurement method based on UAV imagery. The method includes the following steps: extracting two-dimensional farmland boundaries by filtering pixels using frequency domain data and baseline phase parameters; calculating angular deviation mapping to generate a dynamic topological weight coefficient set; introducing an adjustment model with error constraints to optimize and obtain three-dimensional point cloud and orientation data; and extracting profile sequences to perform elevation difference calculations to obtain variation features and filter abrupt change points to output measurement records. In this invention, by combining prior error constraints for adjustment optimization, terrain interference is effectively overcome and systematic ranging offsets of the equipment are eliminated. Elevation gradient calculations are performed on the point cloud, and a penalty term is constructed using roughness and vegetation attributes to deeply suppress high-frequency noise caused by crop shading and gravel. This accurately distinguishes between real embankment fractures and vegetation artifacts, providing an anti-interference judgment benchmark for extracting boundary abrupt change points and ensuring the spatial positioning accuracy of the measurement records.

Description

Technical Field

[0001] This invention relates to the field of photogrammetry, and in particular to a land boundary measurement method based on UAV imagery. Background Technology

[0002] Photogrammetry is a scientific and technological field that extracts geometric and physical information about a target by acquiring photographic images. Its core aspects include image acquisition, image control point deployment, aerial triangulation, and the generation of digital orthophoto maps. Generally, this field utilizes optical sensors mounted on aerial or space platforms to capture two-dimensional images of the earth's surface and features. By analyzing collinearity equations and geometric projection relationships, and combining in-camera and out-of-camera orientation elements, image correction, feature matching, and coordinate transformation are performed to convert the two-dimensional image information into three-dimensional spatial coordinates and extract the spatial position and shape attributes of the target object. Traditional land boundary measurement methods based on UAV imagery refer to the process of using UAVs to acquire aerial images of the survey area and combining them with ground control points to determine land ownership boundaries. The core technical aspect of this method is the determination of the spatial position of land ownership boundary points and boundary lines. In practice, UAVs equipped with visible light cameras typically fly along a preset flight path, taking timed exposures to acquire aerial images with overlap. Simultaneously, ground personnel use real-time dynamic differential equipment to collect the coordinates of control points within the survey area. Subsequently, the aerial images and control point coordinates are input into the photogrammetric processing terminal, and feature point extraction, image matching, and bundle adjustment are performed in sequence to solve the exterior orientation elements of each image. Then, a dense point cloud is generated and a digital surface model is constructed, based on which a digital orthophoto map is output. Finally, the operators manually identify the boundaries of features such as field ridges and ditches on the orthophoto map, and manually draw line segments and closed polygons along the boundaries to extract the boundary point coordinates and generate land boundary lines.

[0003] Traditional boundary measurement relies heavily on manual identification of boundaries and delineation of coordinates on orthophotos. This method is susceptible to interference from the complex environment of farmland. When crops obscure the boundary or gravel is unevenly distributed, it is difficult to identify the actual physical location of the embankment fracture based on the appearance of the image. Vegetation artifacts or surface fragments may be misjudged as ownership boundaries, resulting in random deviations in the spatial location of the boundary. In addition, traditional adjustment calculation processes cannot eliminate the inherent ranging offset of the equipment, which exacerbates the continuous accumulation of positioning errors and makes it difficult for the boundary measurement results to meet the topological connectivity standard. Summary of the Invention

[0004] To address the technical problems existing in the prior art, this invention provides a land boundary measurement method based on UAV imagery.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a land boundary measurement method based on UAV imagery, comprising the following steps: S1: Perform multi-scale frequency domain filtering on farmland surface image sequences, sum the absolute values ​​of the amplitudes of the high-frequency and low-frequency response arrays to obtain high-frequency energy features and low-frequency energy features, and combine them to generate a frequency domain energy dataset; S2: Calculate the ratio of low-frequency to high-frequency energy features in the frequency domain energy dataset to generate texture interference parameters, multiply them with preset base phase parameters to generate dynamic filtering limits, filter farmland surface pixel coordinates whose absolute phase parameters are greater than the dynamic filtering limits and perform line connection processing to generate a two-dimensional farmland boundary set. S3: Calculate the absolute angle deviation parameter between the intersecting line segment parameter and the orthogonal reference angle in the two-dimensional farmland boundary set, sum it with the preset regularization smoothing term, calculate the ratio of the result to the preset information parameter and map and transform it to generate a dynamic topological weight coefficient set; S4: Calculate the product of the line observation error parameter and the dynamic topology weight coefficient set to generate a weighted random parameter, perform spatial forward intersection on the two-dimensional farmland boundary set to obtain the initial three-dimensional positioning parameter, and perform global optimization on the initial three-dimensional positioning parameter in combination with the weighted random parameter to generate three-dimensional point cloud and orientation aggregation data. S5: Along the orthogonal spatial vector corresponding to the two-dimensional boundary parameters of the three-dimensional point cloud and the direction aggregation data, extract the profile sequence from the three-dimensional point cloud parameters, perform first-order differential gradient calculation of elevation to obtain variation feature parameters, filter the coordinates of spatial abrupt change points that are greater than the preset physical sill fracture boundary, and output the land boundary measurement record.

[0006] As a further embodiment of the present invention, the frequency domain energy dataset includes a spectral amplitude matrix, a band energy distribution map, and frequency response extreme points; the two-dimensional farmland boundary set specifically includes an edge pixel cluster, a line segment endpoint coordinate table, and a topological connection diagram; the dynamic topological weight coefficient set specifically includes a confidence scalar vector, a smoothing penalty term matrix, and an angle error compensation value; the three-dimensional point cloud and orientation aggregation data includes spatial coordinate tuples, a normal vector attribute field, and a principal axis orientation vector; and the land boundary measurement record specifically includes a list of ownership boundary nodes, a plot area statistics table, and a sill elevation difference report.

[0007] As a further aspect of the present invention, the steps for obtaining the frequency domain energy dataset are specifically as follows: S101: Acquire the farmland surface image sequence taken by the UAV, perform multi-scale frequency domain filtering on the farmland surface image sequence, extract the corresponding spatial frequency band response matrix, separate the spatial high-frequency signal component and low-frequency signal component, and generate high-frequency response array and low-frequency response array. S102: Based on the high-frequency response array and the low-frequency response array, extract the absolute value of the amplitude for each pixel node of the high-frequency response array and perform matrix element summation to obtain high-frequency energy features; extract the absolute value of the amplitude for each pixel node of the low-frequency response array and perform linear accumulation calculation to obtain low-frequency energy features; and output the energy feature processing result. S103: Call the energy feature processing results, perform data channel splicing operation on the spatial scale dimension for high-frequency energy features and low-frequency energy features, establish a multi-channel energy feature vector, and perform data format encapsulation processing according to multi-scale frequency band attributes to generate a frequency domain energy dataset.

[0008] As a further aspect of the present invention, the steps for obtaining the two-dimensional farmland boundary set are specifically as follows: S201: Call the frequency domain energy dataset, parse the data channels contained in the frequency domain energy dataset, obtain low-frequency energy features and high-frequency energy features, perform division operation to obtain the numerical ratio, establish texture interference parameters, obtain preset basis phase parameters, perform multiplication operation on the basis phase parameters and texture interference parameters to obtain the product, and generate dynamic screening limits. S202: Obtain a farmland surface image sequence, perform phase analysis operation in the two-dimensional frequency domain for each pixel node in the farmland surface image sequence, obtain the absolute phase parameter of the pixel, call the dynamic filtering limit, compare the value of the absolute phase parameter of the pixel with the value of the dynamic filtering limit, filter the pixel nodes whose absolute phase parameter of the pixel is greater than the dynamic filtering limit, extract the planar spatial position information of the filtered pixel nodes, and obtain the coordinates of the farmland surface pixels. S203: Perform planar spatial connectivity analysis on the coordinates of the farmland surface pixels, extract a discrete point set composed of the coordinates of the farmland surface pixels, perform local linear structure analysis and line feature extraction operations on the discrete point set to obtain edge line segment features, perform neighborhood matching and geometric connection processing on the endpoints of the edge line segment features, splice them to form a continuous linear topological structure, and generate a two-dimensional farmland boundary set.

[0009] As a further aspect of the present invention, the process of obtaining the preset substrate phase parameters is specifically as follows: Acquire a set of reference images of farmland without vegetation cover, and perform a two-dimensional frequency domain transformation operation to extract the complex spectrum matrix; Separate the imaginary and real components of the complex spectrum matrix, and perform an arctangent function mapping operation on the imaginary and real components to obtain the basis reference phase map; The phase angle values ​​of pixel nodes in the base reference phase map are extracted and spatial mean calculation is performed to obtain the global phase expectation index; The background thermal noise variance record is extracted by calling the sensor calibration file, and the standard deviation of the background thermal noise is obtained by performing the square root operation and used as the background random phase offset value. The global phase expectation index and the background random phase offset value are added to obtain the initial phase calibration benchmark. The basic orthogonal phase space range is constructed by calling the boundary values ​​of the polar coordinate domain. Interval compression is performed on the initial phase calibration reference term and mapped to the basic orthogonal phase space range to generate the basis phase parameters.

[0010] As a further aspect of the present invention, the step of obtaining the dynamic topology weight coefficient set specifically includes: S301: Perform spatial topological relationship analysis on the linear elements contained in the two-dimensional farmland boundary set, detect the geometric intersection status between line segments, filter the line segment features that are associated with entity intersection, obtain the parameters of the intersecting line segments, perform direction vector inner product and inverse trigonometric function calculation on the parameters of the intersecting line segments, obtain the angle value between the intersecting entities, and generate the plane angle parameter. S302: Based on the plane angle parameter, obtain the preset orthogonal reference angle, perform numerical comparison and subtraction difference operation on the plane angle parameter and the orthogonal reference angle, extract the corresponding angle deviation difference, perform absolute value extraction operation on the angle deviation difference, and generate absolute angle deviation parameter. S303: Call the absolute angle deviation parameter, obtain the preset regularization smoothing term and the preset confidence parameter, perform an addition operation on the absolute angle deviation parameter and the regularization smoothing term to obtain the summation result and use it as the denominator of the division, use the confidence parameter as the numerator of the division to perform the quotient calculation operation, obtain the topology ratio term, perform a parameter mapping transformation operation on the numerical interval of the topology ratio term, and generate a dynamic topology weight coefficient set.

[0011] As a further aspect of the present invention, the process of obtaining the preset orthogonal reference angle is specifically as follows: Call the pre-stored historical standard farmland grid vector data, extract the corresponding boundary interior angle set, count the distribution frequency of the boundary interior angle set, generate a distribution histogram, and set the center angle value corresponding to the peak value of the maximum frequency of the distribution histogram as the preset orthogonal reference angle; The process of obtaining the preset regularization smoothing term is as follows: Read the spatial resolution parameters of the farmland surface image sequence, perform radian conversion mapping between the pixel positioning tolerance record output by the line feature extraction operation and the spatial resolution parameters, obtain the basic angle error bias term and use it as a preset regularization smoothing term; The process of obtaining the preset confidence parameter is as follows: Extract the recognition probability sequence associated with the two-dimensional farmland boundary set, perform mean calculation to obtain the global recognition expectation value, parse the farmland surface image sequence to obtain the signal-to-noise ratio parameter, and multiply the signal-to-noise ratio parameter with the global recognition expectation value to generate a preset confidence parameter.

[0012] As a further aspect of the present invention, the steps for obtaining the three-dimensional point cloud and orientation aggregation data are specifically as follows: S401: Call the dynamic topology weight coefficient set and the two-dimensional farmland boundary set, extract the associated two-dimensional spatial line observation error parameter for the two-dimensional farmland boundary set, perform multiplication calculation of the corresponding elements of the line observation error parameter and the dynamic topology weight coefficient set, combine the error attribute and the topology weight value, and generate a weighted random parameter. S402: Based on the two-dimensional farmland boundary set, extract matching homonymous features of line segment entities from image data of multiple observation perspectives, perform geometric projection intersection operation on spatial beams with matching homonymous features, perform spatial forward intersection processing operation, calculate the initial three-dimensional coordinate values ​​of the spatial entity to be measured, and obtain the initial three-dimensional positioning parameters. S403: Call the weighted random parameters and the initial 3D positioning parameters, construct a diagonal prior weight matrix based on the weighted random parameters, and introduce the prior weight matrix as an error distribution constraint into the bundle adjustment operation model. Perform a global nonlinear iterative optimization calculation operation on the initial 3D positioning parameters, correct the spatial coordinate error values, extract the spatial 3D extension direction vector of the boundary entity, and generate 3D point cloud and orientation aggregation data.

[0013] As a further aspect of the present invention, the steps for obtaining the land boundary measurement records are specifically as follows: S501: Call the three-dimensional point cloud and orientation aggregation data, extract the three-dimensional point cloud parameters and the two-dimensional boundary orientation parameters, generate the orthogonal space vector corresponding to the two-dimensional boundary orientation parameters, and extract the profile data nodes in the three-dimensional point cloud parameters along the orthogonal space vector to obtain the profile sequence. S502: Perform first-order difference gradient calculation on the profile sequence to obtain variation feature parameters; S503: Call the variation feature parameter to obtain the preset physical embankment fracture boundary, compare the variation feature parameter with the physical embankment fracture boundary, filter the coordinates of spatial abrupt change points with values ​​greater than the physical embankment fracture boundary, summarize the coordinates of spatial abrupt change points, perform data assembly processing, and output the land boundary measurement record.

[0014] As a further aspect of the present invention, the process of obtaining the preset solid sill fracture boundary specifically includes: The farmland topographic feature sample database is called to extract a reference embankment three-dimensional point cloud set containing known fault annotation records. The first-order difference gradient of elevation is calculated for the profile nodes in the reference embankment three-dimensional point cloud set to obtain the benchmark sample variation parameter sequence. For the baseline sample variation parameter sequence, perform probability density statistical accounting and fit a distribution curve model, and extract the mathematical expectation mean and discrete standard deviation corresponding to the distribution curve model; The statistical confidence lower bound index is obtained by subtracting the mean of the mathematical expectation from the discrete standard deviation. The basic elevation fluctuation tolerance parameter is extracted by reading the historical geological survey records of the farmland area to be measured. The statistical confidence lower bound index and the basic elevation fluctuation tolerance parameter are added together, and the result of the addition is set as the preset physical embankment fracture limit.

[0015] The beneficial effects of the technical solutions provided in the embodiments of the present invention include at least the following: By combining frequency domain filtering and phase analysis to extract pixel features and fusing topological weights to generate weighted random parameters, prior error constraints are introduced into the bundle adjustment model to perform optimization calculations. This effectively overcomes interference from gentle terrain and eliminates systematic ranging offsets from the measuring equipment. For the 3D point cloud profile sequence, first-order differential gradient calculation of elevation is performed. Nonlinear penalty terms are constructed by introducing surface roughness and vegetation cover attributes. In complex surface environments, high-frequency noise caused by crop shading and gravel is deeply suppressed. The model accurately distinguishes between real ground embankment fractures and vegetation artifacts, providing a high anti-interference judgment benchmark for extracting boundary spatial abrupt change points and ensuring the topological connectivity and spatial positioning accuracy of the final output boundary measurement record. Attached Figure Description

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

[0017] Figure 1 This is a schematic diagram of the workflow of the present invention; Figure 2 This is a detailed flowchart of S1 of the present invention; Figure 3 This is a detailed flowchart of the S2 process of the present invention; Figure 4 This is a detailed flowchart of the S3 process of the present invention; Figure 5 This is a detailed flowchart of the S4 process of the present invention; Figure 6 This is a detailed flowchart of S5 of the present invention. Detailed Implementation

[0018] The technical solution of the present invention will now be described with reference to the accompanying drawings.

[0019] In embodiments of the present invention, words such as "exemplarily," "for example," etc., are used to indicate that something is an example, illustration, or description. Any embodiment or design described as "exemplary" in the present invention should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of the word "exemplary" is intended to present the concept in a concrete manner. Furthermore, in embodiments of the present invention, the meaning expressed by "and / or" can be both, or either one.

[0020] In the embodiments of this invention, the terms "image" and "picture" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, they convey the same meaning. Similarly, the terms "of," "corresponding (relevant)," and "corresponding" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, they convey the same meaning.

[0021] In this embodiment of the invention, sometimes a subscript such as W1 may be written in a non-subscript form such as W1. When the difference is not emphasized, the meaning they express is the same.

[0022] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0023] Please see Figure 1 This invention provides a technical solution: a land boundary measurement method based on UAV imagery, comprising the following steps: S1: Perform multi-scale frequency domain filtering on farmland surface image sequences, sum the absolute values ​​of the amplitudes of the high-frequency and low-frequency response arrays to obtain high-frequency energy features and low-frequency energy features, and combine them to generate a frequency domain energy dataset; S2: Calculate the ratio of low-frequency to high-frequency energy features in the frequency domain energy dataset to generate texture interference parameters, multiply them with preset base phase parameters to generate dynamic filtering limits, filter farmland surface pixel coordinates whose absolute phase parameters are greater than the dynamic filtering limits and perform line connection processing to generate a two-dimensional farmland boundary set. S3: Calculate the absolute angular deviation parameter between the intersecting line segment parameters and the orthogonal reference angle in the two-dimensional farmland boundary set, sum it with the preset regularization smoothing term, calculate the ratio of the result to the preset information parameter and map and transform it to generate a dynamic topological weight coefficient set; S4: Calculate the product of the line observation error parameter and the dynamic topology weight coefficient set to generate a weighted random parameter. Perform spatial forward intersection on the two-dimensional farmland boundary set to obtain the initial three-dimensional positioning parameter. Combine the weighted random parameter to perform global optimization on the initial three-dimensional positioning parameter to generate three-dimensional point cloud and orientation aggregation data. S5: Along the orthogonal spatial vector corresponding to the two-dimensional boundary parameters of the three-dimensional point cloud and the direction aggregation data, extract the profile sequence from the three-dimensional point cloud parameters, perform first-order differential gradient calculation of elevation to obtain variation feature parameters, filter the coordinates of spatial abrupt change points that are greater than the preset physical sill fracture boundary, and output the land boundary measurement record.

[0024] The frequency domain energy dataset includes a spectral amplitude matrix, a band energy distribution map, and frequency response extreme points. The two-dimensional farmland boundary set specifically includes edge pixel clusters, line segment endpoint coordinate tables, and topological connection diagrams. The dynamic topological weight coefficient set specifically includes a confidence scalar vector, a smoothing penalty term matrix, and angle error compensation values. The three-dimensional point cloud and orientation aggregation data includes spatial coordinate tuples, normal vector attribute fields, and principal axis orientation vectors. The land boundary measurement records specifically include a list of ownership boundary nodes, a plot area statistics table, and a sill elevation difference report.

[0025] Please see Figure 2 The specific steps for obtaining the frequency domain energy dataset are as follows: S101: Acquire the farmland surface image sequence taken by the UAV, perform multi-scale frequency domain filtering on the farmland surface image sequence, extract the corresponding spatial frequency band response matrix, separate the spatial high-frequency signal component and low-frequency signal component, and generate high-frequency response array and low-frequency response array. A sequence of farmland surface images captured by a drone was acquired. The red, green, and blue color physical channels of each frame were read according to a pre-set overlap of 80% and a relative flight altitude of 120 meters. For the farmland surface image sequence, the original three-channel color pixel data was converted into a single-channel 8-bit grayscale matrix using a weighted average method. The single-frame resolution was uniformly resampled to 4000 rows by 3000 columns. For the grayscale matrix, a two-dimensional fast Fourier transform algorithm was used to map the spatial domain illumination pixels to the frequency domain, outputting a complex spectrum matrix containing real and imaginary parts. Multi-scale frequency domain filtering was performed on the complex spectrum matrix. A Butterworth low-pass filter mask with an order of 2 and a cutoff frequency parameter set to 50 was constructed. This mask was multiplied by the spectrum matrix to filter out high-frequency edge signals and extract the corresponding low-frequency spatial signal components. A high-pass filter mask complementary to the low-pass filter is constructed, with values ​​of 0 in the central region and 1 in the outer edge region. This high-pass mask is then multiplied element-wise with the original spectral matrix to extract the corresponding spatial high-frequency signal components. Two-dimensional inverse Fourier transforms are performed on the separated spatial high-frequency and low-frequency signal components respectively, remapping the frequency domain complex signals back to the original spatial domain. This generates a high-frequency response array containing details of surface contour fractures and a low-frequency response array containing information about the smooth surface background. Both response arrays maintain a strictly 4000x3000 32-bit floating-point data structure.

[0026] S102: Based on the high-frequency response array and the low-frequency response array, extract the absolute value of the amplitude for each pixel node of the high-frequency response array and perform matrix element summation to obtain high-frequency energy features; extract the absolute value of the amplitude for each pixel node of the low-frequency response array and perform linear accumulation calculation to obtain low-frequency energy features; and output the energy feature processing results. Based on the extracted high-frequency and low-frequency response arrays, for each pixel node of the high-frequency response array, the real and imaginary values ​​are extracted according to the complex form of its inverse transform output. The square root of the sum of the squares of the real and imaginary parts is calculated to extract the absolute amplitude of each coordinate point. Centered on the currently processed high-frequency pixel, a 5x5 two-dimensional sliding window is constructed inside the matrix. The 25 neighboring pixel nodes covered by this window are traversed, the absolute values ​​of the amplitudes of all neighboring nodes are extracted, and the matrix elements are summed. The final sum is directly assigned to the central pixel node. The window with a step size of 1 is traversed sequentially to obtain the high-frequency energy features. For each pixel node in the low-frequency response array, the absolute value of its amplitude in the spatial complex state is also extracted. A pre-constructed Gaussian convolution kernel window is used to replace the basic mean window. The size of the Gaussian kernel is set to 5x5 and the standard deviation parameter is set to 1.2. The absolute value of amplitude within the window is linearly accumulated using this Gaussian kernel, giving the center pixel a higher weight ratio, thereby obtaining smoother low-frequency energy features. After the process is completed, the energy feature processing results are output to the computation control flow. Table 1 shows the detailed configuration of the sliding window and convolution parameters in the energy feature extraction process.

[0027] Table 1. Energy Characteristic Processing Parameters; By precisely defining the parameter matrix and performing kernel function accumulation operations, we can effectively suppress local energy anomalies caused by sudden changes in illumination, ensuring that the energy feature processing results accurately quantify the energy distribution of the surface physical texture.

[0028] S103: Call the energy feature processing results, perform data channel splicing operation on the spatial scale dimension for high-frequency energy features and low-frequency energy features, establish a multi-channel energy feature vector, encapsulate the data format according to the multi-scale frequency band attributes, and generate a frequency domain energy dataset. The energy feature processing results generated in memory are retrieved, and a three-dimensional empty tensor data structure with an initial spatial dimension of 4000 rows, 3000 columns, and 2 data channels is initialized in the main storage area. For high-frequency and low-frequency energy features, a data channel splicing operation is performed along the spatial scale dimension. Specifically, the high-frequency energy feature matrix containing information on abrupt changes in terrain edges is written row by row into the 0th data channel of the three-dimensional tensor, while the low-frequency energy feature matrix containing information on the gentle undulations of the terrain base is simultaneously written into the 1st data channel of the three-dimensional tensor, establishing a multi-channel energy feature vector with strict two-dimensional spatial coordinate alignment. To eliminate computational bias caused by excessively large absolute numerical ranges between high-frequency and low-frequency energy, all data nodes of the multi-channel energy feature vector are traversed to find the global maximum and minimum values. For each floating-point value, linear extremum normalization mathematical calculation is performed to uniformly compress the data within all channels to a standardized numerical range of 0.000 to 1.000. Data format encapsulation is performed based on multi-scale frequency band attributes. Customized file metadata tags are appended to the header of the normalized tensor data packet, including a standard identifier with geographic coordinate system specification code 4326 and calibration data with a spatial sampling physical interval of 0.05 meters. Finally, the three-dimensional tensor carrying the header file is serialized into a binary byte stream format to generate a frequency domain energy dataset that can be called across levels.

[0029] Please see Figure 3 The specific steps for obtaining the set of two-dimensional farmland boundaries are as follows: S201: Call the frequency domain energy dataset, parse the data channels contained in the frequency domain energy dataset, obtain low-frequency energy features and high-frequency energy features, perform division operation to obtain the numerical ratio, establish texture interference parameters, obtain the preset basis phase parameters, perform multiplication operation on the basis phase parameters and texture interference parameters to obtain the product, and generate dynamic screening limits. The process involves retrieving the frequency domain energy dataset stored on disk, parsing the data channel mapping indexes in the dataset's header structure, and precisely retrieving the low-frequency energy feature array stored in channel 1 and the high-frequency energy feature array stored in channel 0 by index. Parallel computing threads are used to traverse the image matrix nodes, extracting spatially aligned low-frequency and high-frequency energy features one by one and performing division operations to obtain the numerical ratio. To prevent hardware anomalies caused by division overflow due to high-frequency energy features approaching zero, a very small smoothing parameter with a constant value of 0.001 is forcibly added to the denominator of the division operation. The division outputs of all pixel nodes are filled into a blank two-dimensional array to establish a texture interference parameter. A preset base phase parameter is obtained; its base value is hard-coded to 0.85 based on standard geostatistics, representing the baseline electromagnetic scattering phase shift characteristics of bare soil surfaces under interference-free conditions. For the scalar of the base phase parameter with a value of 0.85, multiplication operations are performed one by one at the corresponding spatial coordinate positions to obtain the product of the scalar and the matrix element. A threshold truncation correction operation is performed on each numerical point in the global product matrix. The upper logical threshold is set to 3.14 and the lower logical threshold is set to 0.00. All product values ​​that exceed the constraint range are forcibly set to the corresponding boundary thresholds, generating a dynamic filtering boundary matrix.

[0030] S202: Acquire farmland surface image sequence, perform phase analysis operation in two-dimensional frequency domain for each pixel node in the farmland surface image sequence, obtain the absolute phase parameter of the pixel, call the dynamic filtering limit, compare the value of the absolute phase parameter of the pixel with the value of the dynamic filtering limit, filter the pixel nodes whose absolute phase parameter of the pixel is greater than the dynamic filtering limit, extract the planar spatial position information of the filtered pixel nodes, and obtain the coordinates of the farmland surface pixels. The system acquires a sequence of raw farmland surface images from an external device, extracts the single-channel grayscale image matrix, and performs phase analysis in the two-dimensional frequency domain for each pixel node in the farmland surface image sequence. The discrete spatial domain illumination intensity matrix is ​​transformed and reconstructed into a complex frequency matrix using the spatial Hilbert transform function. For each frequency domain two-dimensional coordinate point generated by the analysis, the corresponding real and imaginary components are separated. The system calls the underlying arctangent trigonometric function mathematical library instructions to calculate the inverse trigonometric mapping of the ratio of the imaginary to the real components, outputting phase angle data between -3.14 and +3.14. Then, by extracting the absolute value, it rigidly maps this data to a positive value range of 0.00 to 3.14, obtaining the pixel absolute phase parameter containing the physical texture structure details of each independent coordinate pixel. The system calls a pre-resident dynamic filtering boundary matrix in memory, traverses the entire data area, and compares the pixel absolute phase parameter of the current coordinate point with the value of the dynamic filtering boundary at the same coordinate using an arithmetic logic unit. A binary logical image mask with all initial values ​​of 0 is constructed. If the phase parameter of the current coordinate is strictly greater than the corresponding dynamic filtering limit, the value at the corresponding position of the logical mask is forcibly modified to 1, filtering pixel nodes whose absolute phase parameter is greater than the dynamic filtering limit. The mask is then iterated through all valid node elements with a constant value of 1 to extract the planar spatial position information of the filtered pixel nodes. Combined with the spatial resolution parameters of the image configuration file, this information is converted into the real geographical location in the projected coordinate system using a camera perspective affine transformation matrix, thus obtaining a discrete set of farmland surface pixel coordinates.

[0031] S203: Perform planar spatial connectivity analysis on the coordinates of farmland surface pixels, extract a discrete point set composed of farmland surface pixel coordinates, perform local linear structure analysis and line feature extraction operations on the discrete point set, obtain edge line segment features, perform neighborhood matching and geometric connection processing on the endpoints of the edge line segment features, splice them to form a continuous linear topological structure, and generate a two-dimensional farmland boundary set. A planar spatial connectivity analysis was performed on the acquired set of farmland surface pixel coordinates. An eight-neighbor connected component traversal labeling algorithm was introduced to process all discrete coordinate two-dimensional matrices. A minimum connected pixel clustering area threshold of 50 was set, and a non-maximum suppression operation was performed to remove isolated connected regions with an area less than 50 caused by scattered weeds or soil patches, extracting a core discrete point set with significant geometrical continuity. Local linear structure analysis and line feature extraction were performed on this discrete point set. A random sampling consensus algorithm framework was used to repeatedly randomly sample small batches of subsets from the discrete point set to fit a linear model. The maximum orthogonal distance tolerance threshold from a point to the fitted line was set to 3 pixels, and the maximum number of iterations was set to 500. Based on the goodness-of-fit judgment mechanism, the collinear point set that best matches the surface edge contour was retained, and the horizontal and vertical coordinates of the starting and ending points of the line were written into memory to obtain a series of non-intersecting local edge line segment features. Neighborhood matching and geometric connection processing are performed on the endpoints of edge line segments. A circular constrained search neighborhood with a fixed radius of 15 pixels is constructed, centered on each detected independent endpoint. The spatial two-dimensional Euclidean distance between the current endpoint and other endpoints within the neighborhood, as well as the angle between the extensions of their respective line segments, are calculated. If the distance parameter is strictly less than 15 pixels and the angle difference is less than 15 degrees, a complete point matrix is ​​generated between these two undetermined endpoints using a linear interpolation algorithm. This matrix is ​​then stitched together to form a closed linear topological structure with continuous topological orientation. The final output is a set of two-dimensional farmland boundaries containing all contour boundaries.

[0032] Please see Figure 4 The specific steps for obtaining the dynamic topology weight coefficient set are as follows: S301: Perform spatial topological relationship analysis on the linear elements contained in the two-dimensional farmland boundary set, detect the geometric intersection status between line segments, filter the line segment features that are associated with entity intersection, obtain the parameters of the intersecting line segments, perform direction vector inner product and inverse trigonometric function calculation on the parameters of the intersecting line segments, obtain the angle values ​​between the intersecting entities, and generate the plane angle parameter. Spatial topological relationships are analyzed for numerous irregular linear elements contained in a two-dimensional farmland boundary set. A scanline collision detection algorithm is used to construct an event execution queue for sorting the x-coordinates of line segment nodes, and the entire process is then analyzed. Figure 2The geometric intersection state between line segments is determined. All potential collision bounding boxes of line segments are traversed and examined. Features of intersecting line segments with entity intersection are selected, and the exact spatial coordinates of the two-dimensional plane intersection point are recorded. The geometric position matrices of the starting and ending points of the two source line segments at this intersection point are extracted backwards to obtain the core parameters of the intersecting line segments. The direction vector inner product and inverse trigonometric function calculations are performed on the intersecting line segment parameters, transforming the two line segment entities at the intersection point into a first normalized direction vector and a second normalized direction vector diverging from the origin of the intersection point to both ends. It is assumed that the horizontal and vertical coordinate components of the first direction vector are 0.80 and 0.60, respectively, and the horizontal and vertical coordinate components of the second direction vector are -0.60 and 0.80, respectively. The dot product analytical operation is performed on the above two length-normalized vectors, i.e., multiplying the corresponding coordinate components and summing them to solve for their standard cosine value. The processor calls the underlying inverse cosine mathematical function instruction to calculate the absolute radian value corresponding to the cosine value, and then multiplies it by the constant 180 and divides it by the approximate value of pi, 3.14159, to convert it into a degree measurement unit in a continuous range from 0.0 to 180.0 degrees. The angle values ​​between intersecting entities are then obtained, generating a planar angle parameter that reflects the shape of the angle between farmland plots, and pushing it into the underlying array queue in sequence.

[0033] S302: Based on the plane angle parameter, obtain the preset orthogonal reference angle, perform numerical comparison and subtraction difference calculation between the plane angle parameter and the orthogonal reference angle, extract the corresponding angle deviation difference, perform absolute value extraction operation on the angle deviation difference, and generate absolute angle deviation parameter; Based on the planar angle parameters extracted from the array queue, the underlying hard-coded orthogonal reference angle is obtained. This orthogonal reference angle originates from the ideal grid shape parameters of historical high-standard farmland construction specifications, and its absolute value is rigidly set to 90.0 degrees. For each dynamic angle value in the planar angle parameter sequence, a numerical comparison and subtraction difference operation is performed with the fixed orthogonal reference angle. The actual angle value is used as the minuend, and 90.0 degrees is directly subtracted to extract the corresponding angle deviation difference. For each calculated angle deviation difference, an absolute value extraction operation is performed. A mathematical sign stripping instruction forcibly removes the positive and negative polarity attributes of the difference, retaining only the pure amplitude value of its deviation from the orthogonal reference state. Taking the actual planar angle parameter of a farmland inflection point as an example (85.5 degrees), the angle deviation difference obtained after subtracting 90.0 degrees is -4.5 degrees. After the underlying absolute value extraction processing operation, the final output absolute angle deviation parameter value is reset to a positive 4.5 degrees. All processed positive absolute amplitude values ​​will be packaged and combined into a new one-dimensional high-precision floating-point array structure to generate absolute angular deviation parameters covering all geometrically intersecting nodes of the current measurement block, providing a data reference source for subsequent calculation of weights to correct terrain distortion and topological warp.

[0034] S303: Call the absolute angle deviation parameter, obtain the preset regularization smoothing term and the preset confidence parameter, perform an addition operation on the absolute angle deviation parameter and the regularization smoothing term to obtain the summation result and use it as the denominator of the division, use the confidence parameter as the numerator of the division to perform the quotient calculation operation, obtain the topology ratio term, perform parameter mapping transformation operation on the numerical interval of the topology ratio term, and generate a dynamic topology weight coefficient set; The absolute angle deviation parameter in the buffer sequence is retrieved, and a preset regularization smoothing term is obtained. The mathematical constant value of this regularization smoothing term is statically set to 2.0. Simultaneously, a preset confidence parameter is retrieved, which was quantized and calibrated to 0.88 in the previous multispectral image global confidence assessment. A floating-point addition operation is performed on the extracted absolute angle deviation parameter and the regularization smoothing term to obtain the sum, which is then forced as the denominator for the next division operation. For example, for an absolute angle deviation parameter of 4.5 degrees, adding a smoothing factor of 2.0 results in a denominator that is strictly equal to 6.5. The preset confidence parameter of 0.88 is used as the numerator to perform a division quotient calculation, yielding an exact result of approximately 0.135 for 0.88 divided by 6.5, thus obtaining the basic topological ratio term. A parameter mapping transformation operation is performed on the numerical interval of the topological ratio term. The standard logistic sigmoid nonlinear activation function containing the base of the natural constant is used to perform exponential compression calculation on the above linear ratio, strictly mapping and compressing the unbounded ratio term to the probability convergence interval of 0.00 to 1.00, generating a dynamic topological weight coefficient set.

[0035] Table 2. Dynamic Topology Weight Coefficient Mapping Table; Table 2 shows the final topology weight mapping output details for different angle deviations in the core calculation process. This mapping transformation operation assigns topology weights closer to the upper limit threshold to farmland entity nodes with included angles close to the ideal orthogonal 90.0-degree shape.

[0036] Please see Figure 5 The specific steps for obtaining 3D point cloud and orientation aggregation data are as follows: S401: Call the dynamic topology weight coefficient set and the two-dimensional farmland boundary set, extract the associated two-dimensional spatial line observation error parameters for the two-dimensional farmland boundary set, perform multiplication calculation of the corresponding elements for the line observation error parameters and the dynamic topology weight coefficient set, combine the error attributes and the topology weight values, and generate weighted random parameters. The system utilizes a dynamic topological weight coefficient set generated by a memory mapping mechanism and a set of two-dimensional farmland boundaries retained from previous processing. For each discrete line segment entity element in the two-dimensional farmland boundary set, based on the optical distortion intensity of its image block array and the surface sampling resolution corresponding to its flight altitude, the associated two-dimensional spatial line observation error parameter is extracted. The physical baseline value of this line observation error parameter is set to fluctuate linearly with pixel position within the range of 0.02 meters to 0.08 meters. For the current specific boundary, its specific observation error parameter is precisely quantized and assigned a value of 0.05 meters. The extracted line observation error parameter in the memory stack is then multiplied using the standard multiplication of the corresponding feature matrix elements with the dynamic topological weight coefficient set. This involves directly multiplying the aforementioned 0.05-meter error physical baseline by the weight coefficient of 0.533 found in the mapping table for the corresponding node of the line segment. By combining the simple geometric visual error attribute with the topological weight value reflecting local structural stability, a comprehensive correction stochastic index of 0.02665 meters is calculated for this boundary segment. For all hundreds of farmland boundary features extracted in the full-map coordinate system, repeat the above one-to-one numerical scalar multiplication combination process to generate a weighted random parameter one-dimensional floating-point array that is completely isomorphic to the boundary set in the memory structure index dimension. This array is used to quantify and record the reliability parameters of each independent structural segment within the boundary network.

[0037] S402: Based on the set of two-dimensional farmland boundaries, extract matching homonymous features of line segment entities from image data from multiple observation perspectives, perform geometric projection intersection operation on spatial beams with matching homonymous features, perform spatial forward intersection processing operation, calculate the initial three-dimensional coordinate values ​​of the spatial entity to be measured, and obtain the initial three-dimensional positioning parameters. Based on the vertex coordinate information stored in the 2D farmland boundary set, matching homonymous features of boundary line segments are extracted from UAV 2D image data frames with at least three independent observation perspectives where there is effective field-of-view overlap. An epipolar geometric limit constraint algorithm is employed to establish a basic matrix between two images to be matched. Pixel block sliding scan matching is performed along a narrow search band with a set pixel span of 3 pixels. The normalized cross-correlation coefficient index of local 7x7 image blocks is calculated, and high-confidence feature point pairs with cross-correlation values ​​strictly greater than 0.85 are forcibly selected to determine the homonymous pixel array that uniquely represents the edge of farmland features in the multi-view space. Rigorous 3D geometric projection intersection physical operations are performed on the spatial projection beams of all successfully selected matching homonymous features to extract the six-degree-of-freedom external orientation elements of the recording camera in the air, as well as the internal orientation parameters such as the lens focal length and radial distortion coefficient, constructing a complete mathematical equation for multi-view collinearity constraints. Perform spatial forward intersection processing, input all corresponding point parameters into the overdetermined equation system, use the least squares unbiased estimation method to solve the global optimal root of the nonlinear overdetermined equation system in the three-dimensional coordinate system, back-project from the pixel coordinate system of the plane and calculate the initial three-dimensional coordinate values ​​of the physical space entity to be measured, obtain the initial three-dimensional positioning parameters including the geodetic plane coordinates and relative undulation elevation values, and construct the three-dimensional framework of the surveying site.

[0038] S403: Call the weighted random parameters and the initial 3D positioning parameters, construct a diagonal prior weight matrix based on the weighted random parameters, and introduce the prior weight matrix as an error distribution constraint into the bundle adjustment operation model. Perform a global nonlinear iterative optimization solution operation on the initial 3D positioning parameters, correct the spatial coordinate error values, extract the spatial 3D extension direction vector of the boundary entity, and generate 3D point cloud and orientation aggregation data. The system utilizes weighted random parameters generated by a computational library and initial 3D positioning parameters obtained from the aforementioned beam intersection solution. A large diagonal sparse prior weight matrix is ​​constructed based on the reciprocals of the floating-point values ​​in the weighted random parameter array. In this mathematical structure, boundary entities with smaller errors and higher weights are assigned larger prior weight parameters. This prior weight matrix is ​​introduced as a rigid constraint controlling the error propagation distribution into the bundle adjustment matrix operation model, constructing a large nonlinear spatial error equation system with thousands of unknown nodes. A global nonlinear iterative optimization solution is performed on all the initial 3D positioning parameters. The Levinberg-Marquardt hybrid optimization algorithm engine dynamically adjusts the internal damping factor parameter between the steepest gradient descent direction and the Gauss-Newton approximation method. The maximum number of iterations is rigidly constrained to 50, and the error tolerance convergence threshold parameter for iterative solution is precisely locked at 0.001 meters. After multiple rounds of Jacobian matrix differentiation and vector space iterative updates, the 3D spatial coordinate error values ​​are significantly corrected and smoothed. After the coordinate optimization results meet the threshold convergence requirements, the neighborhood 3D point cloud data clusters around all boundary entities are extracted. For each 3D local surface element, the principal component analysis eigenvalue extraction method is applied to calculate its maximum variance principal eigenvector, accurately extracting the 3D extension direction vector of the boundary entity in the surface space, and outputting the 3D point cloud and orientation aggregated data entity.

[0039] Please see Figure 6 The specific steps for obtaining land boundary survey records are as follows: S501: Call the 3D point cloud and orientation aggregation data, extract the 3D point cloud parameters and 2D boundary orientation parameters, generate the orthogonal space vector corresponding to the 2D boundary orientation parameters, and extract the profile data nodes in the 3D point cloud parameters along the orthogonal space vector to obtain the profile sequence. The system retrieves the aggregated 3D point cloud and orientation data from the bundle adjustment global optimization, traverses the core data structure storage tree, and extracts the 3D point cloud parameters and 2D boundary orientation parameters representing the horizontal extension trend of the boundary for each cluster of farmland measurement points. Assuming the 2D boundary orientation parameter of the target boundary segment currently read from system memory is a standardized unit geometric vector, its lateral and longitudinal extension components in the Cartesian coordinate system are recorded as 0.80 and 0.60, respectively. Based on the orthogonality rule of normal projection in analytical geometry, the lateral and longitudinal components of this 2D orientation vector are swapped, and a sign-inverting instruction is forced on one of the swapped components, generating an orthogonal spatial vector that absolutely corresponds to the 2D boundary orientation parameter. Its vector components evolve to -0.60 and -0.80, ensuring that this new vector is strictly perpendicular to the true physical orientation of the current boundary. Using the orthogonal vector in space as the central scanning axis, a depth projection spatial buffer zone is constructed and defined in the 3D digital model. The control processing logic precisely sets the horizontal expansion threshold width of this zone to 2.0 meters. All discrete elevation points and profile data nodes are extracted from the 3D point cloud parameters covered by the orthogonal spatial vector in this buffer zone. The scattered 3D point cloud within the extracted envelope is strictly projected onto the central one-dimensional profile reference vector line. A strict floating-point ascending sorting operation is then performed based on the absolute projection length of each data node from the profile origin endpoint, resulting in a sorted single-line profile sequence.

[0040] S502: Calculate the first-order difference gradient of elevation for the profile sequence, using the following formula: ; Calculate the characteristic parameters of the variation; in, The variational characteristic parameter is obtained by calculating the first-order difference gradient of elevation. The elevation value of the current node is obtained by reading the elevation attribute of the current coordinates from the profile sequence. The elevation value of the preceding node is obtained by reading the elevation attributes of adjacent preceding coordinates in the profile sequence. This is the instrument deviation, obtained by reading the preset sensor factory calibration error record. The node spacing is calculated by determining the Euclidean distance between the current node and its predecessor in the spatial plane. The surface roughness factor is obtained by extracting the degree of discretization of the local normal vectors of the surface point cloud. The vegetation disturbance coefficient is obtained by mapping and transforming multispectral image data of farmland areas. For the one-dimensional profile sequence extracted after distance sorting, the first-order difference gradient of elevation is calculated sequentially for each adjacent node from top to bottom. The obtained elevation values ​​of the current node (105.50), the previous node (105.10), the instrument deviation (0.05), the node spacing (1.00), the surface roughness factor (0.12), and the vegetation disturbance coefficient (0.24) are substituted into the formula. The numerical simulation is then executed in the central processing unit. The process first calculates the absolute difference between the current node elevation (105.50) and the previous node elevation (105.10), yielding a preliminary elevation difference of 0.40. This absolute difference is then subtracted from the factory-specified instrument ranging inherent physical deviation of 0.05, resulting in a net elevation physical change of 0.35. A maximum value discriminator compares 0.35 with the constant 0, and after confirming it is not negative, 0.35 is retained. This net change is then divided by the planar geometric distance between the two nodes (1.00), calculating the basic spatial elevation gradient value as 0.35. Subsequently, the process adds the surface roughness factor (0.12) and the vegetation interference coefficient (0.24) obtained from the extracted spectrum, resulting in a sum of 0.36. A square root operation is then called on 0.36 to calculate the environmental adjustment physical weight factor as 0.60. Multiplying the base spatial elevation gradient of 0.35 by the adjustment weight of 0.60 using a scalar yields a final variation feature parameter of 0.21. This result indicates that the spatial node at the current detection location generates a terrain step phenomenon in the 3D scene that breaks through the conventional background noise. The advantage of this core formula lies in its innovative integration of a nonlinear joint square root penalty weight mechanism reflecting both surface physical properties and biological cover properties into the derivation of the micro-gradient. This mechanism filters out pseudo-abrupt noise points in complex agricultural surface environments with weed interference and gravel tillage characteristics, while preserving and enhancing data representing true fault and ridge features.

[0041] S503: Call the variation feature parameter to obtain the preset solid embankment fracture boundary, compare the variation feature parameter with the solid embankment fracture boundary, filter the coordinates of spatial abrupt change points with values ​​greater than the solid embankment fracture boundary, summarize the coordinates of spatial abrupt change points, perform data assembly processing, and output the land boundary measurement record. The variable feature parameter cache column output by the computational logic is invoked. The preset physical sill fracture boundary is obtained from the pre-mounted business geological feature parameter library configuration file. This boundary physical threshold is derived from a large amount of prior on-site manual profile mapping statistical analysis and is strictly locked at 0.15 in this execution parameter. In the main control register, a comparison instruction is used to compare the variable feature parameter of each detection profile sequence node with the fixed value of the physical sill fracture boundary. A strict greater-than-sign Boolean logic judgment process is executed to eliminate all smooth transition nodes and weak undulation noise points with variable feature parameters lower than or equal to 0.15. The coordinates of deterministic spatial abrupt change points with variable feature parameter values ​​greater than the physical sill fracture boundary are filtered. Since the calculated variable parameter of 0.21 is absolutely greater than the boundary of 0.15, the coordinate location is clearly determined to be an undamaged real physical sill fracture spatial occurrence zone. The horizontal and vertical coordinates and elevation data of the nodes of this occurrence zone in the absolute geodetic control network are extracted. The system aggregates the coordinates of tens of thousands of spatial abrupt change points that meet the threshold conditions from the extraction, analysis, and filtering of the full-domain imagery. It then centrally executes the spatial mapping data assembly and processing mechanism in the memory pool, merging and packaging these data into an object set containing timestamp markers, endpoints of the land rights-confirmed polygons, and standard labels for attribute fields. Finally, it outputs high-precision land boundary measurement records for business terminals.

[0042] Table 3. Encapsulation Structure of Land Boundary Measurement Records; Table 3 shows the detailed physical structure of the underlying fields of the land boundary survey record output to the survey database. The three-dimensional boundary confirmation survey file constructed based on this table can be directly imported into the geographic database of the agricultural asset management database platform, realizing automated end-to-end data flow from contactless surveying to legal property rights delineation.

[0043] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of protection of the described technical solutions.

Claims

1. A land boundary measurement method based on UAV imagery, characterized in that, The method includes: S1: Perform multi-scale frequency domain filtering on farmland surface image sequences, sum the absolute values ​​of the amplitudes of the high-frequency and low-frequency response arrays to obtain high-frequency energy features and low-frequency energy features, and combine them to generate a frequency domain energy dataset; S2: Calculate the ratio of low-frequency to high-frequency energy features in the frequency domain energy dataset to generate texture interference parameters, multiply them with preset base phase parameters to generate dynamic filtering limits, filter farmland surface pixel coordinates whose absolute phase parameters are greater than the dynamic filtering limits and perform line connection processing to generate a two-dimensional farmland boundary set. S3: Calculate the absolute angle deviation parameter between the intersecting line segment parameter and the orthogonal reference angle in the two-dimensional farmland boundary set, sum it with the preset regularization smoothing term, calculate the ratio of the result to the preset information parameter and map and transform it to generate a dynamic topological weight coefficient set; S4: Calculate the product of the line observation error parameter and the dynamic topology weight coefficient set to generate a weighted random parameter, perform spatial forward intersection on the two-dimensional farmland boundary set to obtain the initial three-dimensional positioning parameter, and perform global optimization on the initial three-dimensional positioning parameter in combination with the weighted random parameter to generate three-dimensional point cloud and orientation aggregation data. S5: Along the orthogonal spatial vector corresponding to the two-dimensional boundary parameters of the three-dimensional point cloud and the direction aggregation data, extract the profile sequence from the three-dimensional point cloud parameters, perform first-order differential gradient calculation of elevation to obtain variation feature parameters, filter the coordinates of spatial abrupt change points that are greater than the preset physical sill fracture boundary, and output the land boundary measurement record.

2. The land boundary measurement method based on UAV imagery according to claim 1, characterized in that, The frequency domain energy dataset includes a spectral amplitude matrix, a band energy distribution map, and frequency response extreme points. The two-dimensional farmland boundary set specifically includes edge pixel clusters, a line segment endpoint coordinate table, and a topological connection diagram. The dynamic topological weight coefficient set specifically includes a confidence scalar vector, a smoothing penalty term matrix, and an angle error compensation value. The three-dimensional point cloud and orientation aggregation data includes spatial coordinate tuples, a normal vector attribute field, and a principal axis orientation vector. The land boundary measurement record specifically includes a list of ownership boundary nodes, a plot area statistics table, and a sill elevation difference report.

3. The land boundary measurement method based on UAV imagery according to claim 1, characterized in that, The specific steps for obtaining the frequency domain energy dataset are as follows: S101: Acquire the farmland surface image sequence taken by the UAV, perform multi-scale frequency domain filtering on the farmland surface image sequence, extract the corresponding spatial frequency band response matrix, separate the spatial high-frequency signal component and low-frequency signal component, and generate high-frequency response array and low-frequency response array. S102: Based on the high-frequency response array and the low-frequency response array, extract the absolute value of the amplitude for each pixel node of the high-frequency response array and perform matrix element summation to obtain high-frequency energy features; extract the absolute value of the amplitude for each pixel node of the low-frequency response array and perform linear accumulation calculation to obtain low-frequency energy features; and output the energy feature processing result. S103: Call the energy feature processing results, perform data channel splicing operation on the spatial scale dimension for high-frequency energy features and low-frequency energy features, establish a multi-channel energy feature vector, and perform data format encapsulation processing according to multi-scale frequency band attributes to generate a frequency domain energy dataset.

4. The land boundary measurement method based on UAV imagery according to claim 3, characterized in that, The specific steps for obtaining the set of two-dimensional farmland boundaries are as follows: S201: Call the frequency domain energy dataset, parse the data channels contained in the frequency domain energy dataset, obtain low-frequency energy features and high-frequency energy features, perform division operation to obtain the numerical ratio, establish texture interference parameters, obtain preset basis phase parameters, perform multiplication operation on the basis phase parameters and texture interference parameters to obtain the product, and generate dynamic screening limits. S202: Obtain a farmland surface image sequence, perform phase analysis operation in the two-dimensional frequency domain for each pixel node in the farmland surface image sequence, obtain the absolute phase parameter of the pixel, call the dynamic filtering limit, compare the value of the absolute phase parameter of the pixel with the value of the dynamic filtering limit, filter the pixel nodes whose absolute phase parameter of the pixel is greater than the dynamic filtering limit, extract the planar spatial position information of the filtered pixel nodes, and obtain the coordinates of the farmland surface pixels. S203: Perform planar spatial connectivity analysis on the coordinates of the farmland surface pixels, extract a discrete point set composed of the coordinates of the farmland surface pixels, perform local linear structure analysis and line feature extraction operations on the discrete point set to obtain edge line segment features, perform neighborhood matching and geometric connection processing on the endpoints of the edge line segment features, splice them to form a continuous linear topological structure, and generate a two-dimensional farmland boundary set.

5. The land boundary measurement method based on UAV imagery according to claim 4, characterized in that, The process of obtaining the preset base phase parameters is as follows: Acquire a set of reference images of farmland without vegetation cover, and perform a two-dimensional frequency domain transformation operation to extract the complex spectrum matrix; Separate the imaginary and real components of the complex spectrum matrix, and perform an arctangent function mapping operation on the imaginary and real components to obtain the basis reference phase map; The phase angle values ​​of pixel nodes in the base reference phase map are extracted and spatial mean calculation is performed to obtain the global phase expectation index; The background thermal noise variance record is extracted by calling the sensor calibration file, and the standard deviation of the background thermal noise is obtained by performing the square root operation and used as the background random phase offset value. The global phase expectation index and the background random phase offset value are added to obtain the initial phase calibration benchmark. The basic orthogonal phase space range is constructed by calling the boundary values ​​of the polar coordinate domain. Interval compression is performed on the initial phase calibration reference term and mapped to the basic orthogonal phase space range to generate the basis phase parameters.

6. The land boundary measurement method based on UAV imagery according to claim 4, characterized in that, The specific steps for obtaining the dynamic topology weight coefficient set are as follows: S301: Perform spatial topological relationship analysis on the linear elements contained in the two-dimensional farmland boundary set, detect the geometric intersection status between line segments, filter the line segment features that are associated with entity intersection, obtain the parameters of the intersecting line segments, perform direction vector inner product and inverse trigonometric function calculation on the parameters of the intersecting line segments, obtain the angle value between the intersecting entities, and generate the plane angle parameter. S302: Based on the plane angle parameter, obtain the preset orthogonal reference angle, perform numerical comparison and subtraction difference operation on the plane angle parameter and the orthogonal reference angle, extract the corresponding angle deviation difference, perform absolute value extraction operation on the angle deviation difference, and generate absolute angle deviation parameter. S303: Call the absolute angle deviation parameter, obtain the preset regularization smoothing term and the preset confidence parameter, perform an addition operation on the absolute angle deviation parameter and the regularization smoothing term to obtain the summation result and use it as the denominator of the division, use the confidence parameter as the numerator of the division to perform the quotient calculation operation, obtain the topology ratio term, perform a parameter mapping transformation operation on the numerical interval of the topology ratio term, and generate a dynamic topology weight coefficient set.

7. The land boundary measurement method based on UAV imagery according to claim 6, characterized in that, The process of obtaining the preset orthogonal reference angle is as follows: Call the pre-stored historical standard farmland grid vector data, extract the corresponding boundary interior angle set, count the distribution frequency of the boundary interior angle set, generate a distribution histogram, and set the center angle value corresponding to the peak value of the maximum frequency of the distribution histogram as the preset orthogonal reference angle; The process of obtaining the preset regularization smoothing term is as follows: Read the spatial resolution parameters of the farmland surface image sequence, perform radian conversion mapping between the pixel positioning tolerance record output by the line feature extraction operation and the spatial resolution parameters, obtain the basic angle error bias term and use it as a preset regularization smoothing term; The process of obtaining the preset confidence parameter is as follows: Extract the recognition probability sequence associated with the two-dimensional farmland boundary set, perform mean calculation to obtain the global recognition expectation value, parse the farmland surface image sequence to obtain the signal-to-noise ratio parameter, and multiply the signal-to-noise ratio parameter with the global recognition expectation value to generate a preset confidence parameter.

8. The land boundary measurement method based on UAV imagery according to claim 6, characterized in that, The specific steps for obtaining the 3D point cloud and orientation aggregation data are as follows: S401: Call the dynamic topology weight coefficient set and the two-dimensional farmland boundary set, extract the associated two-dimensional spatial line observation error parameter for the two-dimensional farmland boundary set, perform multiplication calculation of the corresponding elements of the line observation error parameter and the dynamic topology weight coefficient set, combine the error attribute and the topology weight value, and generate a weighted random parameter. S402: Based on the two-dimensional farmland boundary set, extract matching homonymous features of line segment entities from image data of multiple observation perspectives, perform geometric projection intersection operation on spatial beams with matching homonymous features, perform spatial forward intersection processing operation, calculate the initial three-dimensional coordinate values ​​of the spatial entity to be measured, and obtain the initial three-dimensional positioning parameters. S403: Call the weighted random parameters and the initial 3D positioning parameters, construct a diagonal prior weight matrix based on the weighted random parameters, and introduce the prior weight matrix as an error distribution constraint into the bundle adjustment operation model. Perform a global nonlinear iterative optimization calculation operation on the initial 3D positioning parameters, correct the spatial coordinate error values, extract the spatial 3D extension direction vector of the boundary entity, and generate 3D point cloud and orientation aggregation data.

9. The land boundary measurement method based on UAV imagery according to claim 8, characterized in that, The specific steps for obtaining the land boundary measurement records are as follows: S501: Call the three-dimensional point cloud and orientation aggregation data, extract the three-dimensional point cloud parameters and the two-dimensional boundary orientation parameters, generate the orthogonal space vector corresponding to the two-dimensional boundary orientation parameters, and extract the profile data nodes in the three-dimensional point cloud parameters along the orthogonal space vector to obtain the profile sequence. S502: Perform first-order difference gradient calculation on the profile sequence to obtain variation feature parameters; S503: Call the variation feature parameter to obtain the preset physical embankment fracture boundary, compare the variation feature parameter with the physical embankment fracture boundary, filter the coordinates of spatial abrupt change points with values ​​greater than the physical embankment fracture boundary, summarize the coordinates of spatial abrupt change points, perform data assembly processing, and output the land boundary measurement record.

10. The land boundary measurement method based on UAV imagery according to claim 9, characterized in that, The process of obtaining the preset physical sill fracture boundary is as follows: The farmland topographic feature sample database is called to extract a reference embankment three-dimensional point cloud set containing known fault annotation records. The first-order difference gradient of elevation is calculated for the profile nodes in the reference embankment three-dimensional point cloud set to obtain the benchmark sample variation parameter sequence. For the baseline sample variation parameter sequence, perform probability density statistical accounting and fit a distribution curve model, and extract the mathematical expectation mean and discrete standard deviation corresponding to the distribution curve model; The statistical confidence lower bound index is obtained by subtracting the mean of the mathematical expectation from the discrete standard deviation. The basic elevation fluctuation tolerance parameter is extracted by reading the historical geological survey records of the farmland area to be measured. The statistical confidence lower bound index and the basic elevation fluctuation tolerance parameter are added together, and the result of the addition is set as the preset physical embankment fracture limit.

Images