Land area information acquisition system based on territorial space planning

Abstract

The application relates to the technical field of land area measurement, and discloses a land area information acquisition system based on national land space planning, which comprises a data acquisition module, a digital elevation model, a terrace three-coupling perimeter gain calculation module, a direction sensitivity coefficient calculation module, an equivalent outward expansion radius calculation module, an area correction module and a summary output module. The terrace frequency band is extracted according to the digital elevation model, and a terrace weight is generated, so that non-terrace area terrain noise can be filtered. The three-coupling perimeter gain kernel is constructed by fusing the slope projection scaling light visibility and the sensitive items of the orthographic and elevation model errors, so as to meet the high-precision and high-credibility requirements of the national land space planning on the area information of the terrace area.

Description

Technical Field

[0001] This invention relates to the field of land area measurement technology, and more specifically, to a land area information collection system based on national land spatial planning. Background Technology

[0002] Territorial spatial planning places stringent demands on the accuracy and reliability of land area data. Terraced fields, a common land use type in mountainous areas of my country, present unique challenges to area data collection due to their strip-shaped landform characteristics. Traditional land area collection techniques often employ a comprehensive processing approach, first directly analyzing the digital elevation model as a whole, and then combining it with satellite orthophotos to extract plot boundaries and calculate the area. However, this approach fails to specifically extract the geomorphic frequency bands corresponding to the terrace strips, easily introducing topographic noise from non-terraced areas into subsequent measurement processes, interfering with the accurate perimeter and area calculations.

[0003] In the perimeter correction stage, traditional techniques typically only consider the influence of map projection or terrain slope, ignoring the coupled effects of multiple factors such as slope projection scaling, visibility changes caused by lighting and shadows, and boundary jitter caused by vertical errors in digital elevation models under off-axis imaging. This makes it impossible to accurately reflect the actual deviation of the perimeter of terraced plots. Furthermore, UAV sampling often uses uniform grid flight paths, failing to set reasonable sampling intervals based on the wavelength of the terrace strips, and neglecting to target highly sensitive areas for perimeter deviation for intensive sampling. This makes it difficult to obtain high-precision data for physical calibration of key correction parameters.

[0004] Furthermore, traditional techniques lack a clear geometric mapping logic from perimeter gain to area correction and have not established a complete uncertainty propagation system. As a result, the area correction results lack reliable theoretical support, and the accuracy cannot be effectively traced. Ultimately, the land area data of the terraced areas produced has a large deviation, which is difficult to meet the decision-making and control needs of national land spatial planning. Summary of the Invention

[0005] This invention provides a land area information collection system based on national land spatial planning, which solves the technical problems mentioned in the background.

[0006] This invention provides a land area information collection system based on national land spatial planning, comprising:

[0007] The data acquisition module is used to acquire digital elevation models, calculate slope, slope direction, contour line normal, generate terrace weights and terrace frequency band wavelengths through bandpass filtering, and calculate terrain elevation angle and visibility by combining solar geometric parameters.

[0008] The terraced field three-coupling perimeter gain calculation module is used to segment land types in satellite orthophotos to obtain boundary curves, generate length scaling by combining slope, slope direction and projection scale factor, integrate visibility and orthophoto and elevation model error sensitive terms, construct a three-coupling perimeter gain kernel, and obtain the terraced field three-coupling perimeter gain by weighted accumulation of terraced field weights.

[0009] The directional sensitivity coefficient calculation module is used to set the UAV flight path spacing according to the wavelength of the terrace frequency band, determine the key acquisition zone according to the product of the terrace weight and the three-coupled perimeter gain kernel, calculate the directional sensitivity coefficient through UAV data, and update the error sensitivity terms of the orthophoto and elevation models and the three-coupled perimeter gain of the terrace based on this.

[0010] The equivalent outer radius calculation module is used to obtain the frequency band perimeter of the terraced fields through weighted accumulation of terraced field weights, determine the equivalent outer radius based on the ratio of the terraced field three-coupling perimeter gain to the perimeter, and calculate the brightness gradient based on satellite orthophoto and determine the sign based on the solar azimuth angle.

[0011] The area correction module is used to take the product of the sign, the equivalent outer radius, and the perimeter of the terrace frequency band as the area correction amount, subtract the correction amount from the original area to obtain the corrected area, and at the same time calculate the area uncertainty.

[0012] The summary output module is used to summarize the corrected area and corresponding uncertainty of each land type as the final output.

[0013] The beneficial effects of this invention are as follows: This invention extracts terrace frequency bands and generates terrace weights based on a digital elevation model, which can filter terrain noise in non-terraced areas and achieve directional data processing; by fusing slope projection scaling, illumination visibility, and orthophoto and elevation model error sensitivity terms to construct a three-coupled perimeter gain kernel, it accurately quantifies the superimposed influence of multiple physical factors on the perimeter; by combining the wavelength of the terrace frequency band to design an adaptive flight path for UAVs and delineating key acquisition zones, and by utilizing the UAV data calibration direction sensitivity coefficient, the reliability of perimeter gain calculation is further improved. Simultaneously, a clear geometric mapping from perimeter gain to area correction is established, along with a complete uncertainty propagation system, ensuring traceable data accuracy. The final output land use summary data can meet the high-precision and high-reliability requirements of land spatial planning for terraced area area information. Attached Figure Description

[0014] Figure 1 This is a schematic diagram of a land area information collection system based on national land spatial planning according to the present invention.

[0015] In the diagram: Data acquisition module 101, terrace three-coupling perimeter gain calculation module 102, directional sensitivity coefficient calculation module 103, equivalent outward radius calculation module 104, area correction module 105, summary output module 106. Detailed Implementation

[0016] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0017] It should be noted that, unless otherwise defined, the technical or scientific terms used in one or more embodiments of the present invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in one or more embodiments of the present invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" indicate that the element or object preceding the term encompasses the elements or objects listed following the term and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0018] like Figure 1 As shown, a land area information collection system based on territorial spatial planning includes:

[0019] Data acquisition module 101 is used to acquire digital elevation model, calculate slope, slope direction, contour line normal, generate terrace weight and terrace frequency band wavelength through bandpass filtering, and calculate terrain elevation angle and visibility by combining solar geometric parameters;

[0020] The terraced field three-coupling perimeter gain calculation module 102 is used to perform land use segmentation on satellite orthophotos to obtain boundary curves, combine slope, slope direction and projection scale factor to generate length scaling, fuse visibility and orthophoto and elevation model error sensitive terms, construct a three-coupling perimeter gain kernel, and obtain the terraced field three-coupling perimeter gain by weighted accumulation of terraced field weights.

[0021] The directional sensitivity coefficient calculation module 103 is used to set the UAV flight path spacing according to the wavelength of the terrace frequency band, determine the key acquisition zone according to the product of the terrace weight and the three-coupled perimeter gain kernel, calculate the directional sensitivity coefficient through UAV data, and update the error sensitivity terms of the orthophoto and elevation models and the three-coupled perimeter gain of the terrace accordingly.

[0022] The equivalent outer radius calculation module 104 is used to obtain the perimeter of the terrace frequency band by weighted accumulation of terrace weights, determine the equivalent outer radius according to the ratio of the terrace three-coupling perimeter gain to the perimeter, and calculate the brightness gradient according to the satellite orthophoto and determine the sign by combining the solar azimuth angle.

[0023] The area correction module 105 is used to take the product of the sign, the equivalent outer radius, and the perimeter of the terrace frequency band as the area correction amount, subtract the correction amount from the original area to obtain the corrected area, and at the same time calculate the area uncertainty.

[0024] The summary output module 106 is used to summarize the corrected area and corresponding uncertainty of each land type as the final output.

[0025] In one embodiment of the present invention, the first-order partial derivatives of the digital elevation model in the horizontal and vertical directions of the plane coordinate system are calculated, the slope is determined based on the arctangent of the square root of the sum of the squares of the partial derivatives in these two directions, the direction of maximum slope is determined based on the arctangent of the ratio of the partial derivatives in these two directions, and the angle value of the direction of maximum slope is converted into a unit vector of the cross-contour line normal in the two-dimensional plane.

[0026] Centered on each pixel, a local elevation profile is extracted along the unit vector normal to the contour line. The local elevation profile is then transformed in the frequency domain to obtain the local energy spectrum. The frequency corresponding to the maximum energy in the local energy spectrum is identified as the main peak frequency. The reciprocal of the main peak frequency is calculated to obtain the main wavelength of the terrace band. At the same time, the energy spectrum within a preset neighborhood of the main peak frequency is integrated to obtain the bandpass signal energy. The bandpass signal energy is divided by the sum of the integrals of the bandpass signal energy of the entire map to obtain the normalized terrace weight.

[0027] A unit vector for the sun's direction is constructed based on the trigonometric function value of the sun's azimuth angle. The maximum elevation angle of each pixel relative to the surrounding terrain is searched and calculated along this vector to determine the terrain elevation angle. The difference between the sun's elevation angle and the terrain elevation angle is calculated, and this difference is mapped to a value between zero and one using a non-linear smoothing mapping function to obtain the visibility.

[0028] It should be noted that the first-order partial derivatives of the digital elevation model (DEM) in the horizontal and vertical directions of the plane coordinate system represent the rates of elevation change of the DEM in the horizontal and vertical directions, respectively, reflecting the degree of elevation inclination in the corresponding directions. Slope represents the angle of inclination at a point on the Earth's surface, reflecting the steepness of the terrain at that point. The direction of maximum gradient represents the horizontal direction in which the elevation of a point on the surface decreases the fastest, i.e., the location of the steepest downhill slope. The unit vector across contour lines represents a unit-length two-dimensional vector perpendicular to the contour lines and pointing in the direction of maximum gradient, used to define the interception direction of the terrain profile. A local elevation profile represents a continuous sequence of elevation data intercepted along the direction of the unit vector across contour lines, centered on a pixel in the DEM. The local energy spectrum represents the energy distribution corresponding to each frequency component obtained after frequency domain transformation of the local elevation profile (using Fast Fourier Transform), reflecting the spatial frequency characteristics of the profile.

[0029] It should be noted that the main peak frequency represents the frequency component with the highest energy value in the local energy spectrum, corresponding to the most significant spatial rhythm of the terraced stripes. The dominant wavelength of the terraced band represents the spatial wavelength corresponding to the main peak frequency, reflecting the typical width of the terraced stripes. The bandpass signal energy represents the energy integral value within a preset neighborhood of the main peak frequency in the local energy spectrum, reflecting the energy intensity of the terraced band. The integral sum of the bandpass signal energy of the entire image represents the accumulated value of the bandpass signal energy of all pixels in the entire processing area, used for normalizing the terraced weights. The normalized terraced weight represents the ratio of the bandpass signal energy of a single pixel to the integral sum of the bandpass signal energy of the entire image, reflecting the degree to which that point conforms to the characteristics of the terraced stripes.

[0030] It should be noted that the solar azimuth angle represents the angle between the projection of sunlight onto the horizontal plane and true north, reflecting the horizontal orientation of the sun. The solar direction unit vector represents a unit-length vector on the horizontal plane constructed based on the solar azimuth angle, pointing horizontally towards the sun. The terrain elevation angle represents the maximum elevation angle of a point on the Earth's surface looking towards the surrounding terrain along the solar direction unit vector, reflecting the degree to which that point is obscured by the surrounding terrain. The solar elevation angle represents the angle between sunlight and the horizontal plane, reflecting the vertical altitude of the sun. The difference between the solar elevation angle and the terrain elevation angle is used to determine the state of light occlusion at a point on the Earth's surface. The nonlinear smoothing mapping function (Sigmoid) represents a smooth nonlinear function that maps any value to the interval between zero and one, used to soften the binary determination of light occlusion. Visibility represents the value obtained after processing by the nonlinear smoothing mapping function, reflecting the degree of light visibility at a point on the Earth's surface, with a value ranging from zero to one.

[0031] Specifically, the formula for calculating the unit vector across contour lines is to take the cosine of the angle of the direction of maximum slope as the horizontal axis component and the sine as the vertical axis component, and then normalize the two-dimensional vector to a unit length; the formula for calculating the local energy spectrum is to perform a fast Fourier transform on the local elevation profile and take the square of the magnitude of the transform result; the formula for calculating the normalized terrace weight is to take the bandpass signal energy of a single pixel and divide it by the integral sum of the bandpass signal energy of all pixels in the entire image; the formula for calculating the unit vector in the solar direction is to take the cosine of the solar azimuth angle as the horizontal axis component and the sine as the vertical axis component, and then normalize the two-dimensional vector to a unit length; the formula for calculating the terrain elevation angle is to traverse terrain points at different distances along the direction of the solar direction unit vector, calculate the elevation angle of each point relative to the target pixel, and take the maximum value.

[0032] It should be noted that the local elevation profile cut-off length can be centered on the target pixel, extending 5 pixels to each side along the unit vector direction of the contour line normal, resulting in an overall profile length corresponding to 11 pixels of actual ground distance. The preset neighborhood range of the main peak frequency can be set to the frequency range of 90% to 110% of the main peak frequency. The search distance along the unit vector direction of the sun can be set to terrain points within a 30-meter radius around the target pixel.

[0033] In one embodiment of the present invention, the boundary curve of the land parcel is obtained by segmenting and extracting the satellite orthophoto image. The differential of each point on the boundary curve in the vertical axis direction and the differential in the horizontal axis direction are calculated respectively. The arctangent value of the ratio of the two in the four quadrants is taken as the boundary tangential orientation.

[0034] Calculate the square of the cosine of the difference between the boundary tangential azimuth and the maximum slope direction, multiply it by the square of the slope tangent, add one, and take the reciprocal of the square root to obtain the slope projection scaling component. Multiply the slope projection scaling component by the map projection scale factor to obtain the length scaling amount.

[0035] Specifically, the slope projection scaling component at point x on the ground surface. The calculation formula is as follows:

[0036]

[0037] in , and These represent the slope, boundary tangential orientation, and direction of maximum gradient at point x on the ground surface, respectively.

[0038] It should be noted that satellite orthophotos represent surface images projected onto a plane of equal scale after eliminating geometric distortions such as sensor viewpoint and topographic relief. They can be used to accurately extract land parcel boundaries and land use information. The boundary curves of land parcels represent the outlines of different land use parcels after the satellite orthophotos have been segmented by land use. The differentials of each point on the boundary curve along the vertical and horizontal axes represent the rate of change of the vertical and horizontal coordinates of the boundary curve at a given point with respect to arc length, reflecting the slope of the tangent lines in the vertical and horizontal directions at that point. The boundary tangential azimuth represents the horizontal azimuth angle of the tangent line at a given point on the boundary curve, reflecting the direction of the boundary's extension at that point. The square of the cosine of the difference between the boundary tangential azimuth and the direction of maximum slope gradient is used to quantify the parallelism between the boundary and contour lines. The slope projection scaling component represents the linear scaling ratio of the actual slope boundary length projected onto the horizontal plane, reflecting the influence of terrain slope on the boundary length. The map projection scale factor represents the fine-tuning coefficient for the linear scale of the surface during map projection, used to correct scale deviations introduced by the projection. The length scaling factor represents the combined length correction factor that integrates the slope projection scaling and the map projection scale, and can calibrate the planar observation length of the boundary.

[0039] In one embodiment of the present invention, the absolute value of the cosine of the difference between the boundary tangential azimuth and the sensor line-of-sight azimuth is calculated, and multiplied by the direction sensitivity coefficient, the vertical mean square error of the digital elevation model, and the tangent of the off-axis angle. The result of the multiplication is divided by the ground sampling distance, and one is added to the quotient to obtain the orthophoto and elevation model error sensitivity terms.

[0040] Multiply the length scaling, visibility, and orthorectification and elevation model error sensitivity terms to obtain the three-coupled perimeter gain kernel;

[0041] Integrate the product of the three-coupled perimeter gain kernel and the terrace weight along the boundary curve, and simultaneously integrate the terrace weight along the boundary curve. Divide the result of the former by the result of the latter, and subtract one from the quotient to obtain the three-coupled perimeter gain of the terrace.

[0042] Specifically, the orthorectification and elevation model error sensitivity terms at point x on the ground surface. The calculation formula is as follows:

[0043]

[0044] in Indicates the directional sensitivity coefficient. This represents the vertical mean square error of the digital elevation model at point x on the Earth's surface. Indicates the off-axis angle. Indicates the ground sampling distance. Indicates the tangential orientation of the boundary at point x on the Earth's surface. Indicates the sensor's line of sight orientation.

[0045] It should be noted that the sensor line-of-sight azimuth represents the azimuth angle of the projection of the satellite sensor's line of sight onto the horizontal plane when observing the Earth's surface, reflecting the sensor's horizontal observation orientation. The absolute value of the cosine of the difference between the boundary tangential azimuth and the sensor line-of-sight azimuth is used to quantify the sensitivity of the boundary orientation to observation errors. The direction sensitivity coefficient reflects the sensitivity coefficient of converting the vertical error of the digital elevation model into planar boundary jitter. The vertical mean square error of the digital elevation model represents the statistical error of the digital elevation model's elevation data, reflecting the accuracy level of the elevation data. The off-axis angle is the angle between the satellite sensor's line of sight and the ground normal, reflecting the degree of tilt observation by the sensor. The ground sampling distance represents the actual ground size corresponding to a single pixel in the satellite orthophoto image, reflecting the spatial resolution of the image. The orthophoto and elevation model error sensitivity term represents the boundary length deviation coefficient caused by the combined vertical error of the digital elevation model and the sensor's observation angle, used to correct the perimeter distortion caused by observation errors. The three-coupled perimeter gain kernel represents the comprehensive gain coefficient that fuses the length scaling, illumination visibility, and observation error sensitivity term, reflecting the actual degree of change in the boundary perimeter under the superposition of multiple factors. The product of the three-coupled perimeter gain kernel and the terrace weights represents the weighted gain value that limits the overall perimeter gain to the terrace frequency band, ensuring that only the perimeter deviation of the terrace topography is corrected. The terrace three-coupled perimeter gain represents the net gain or net attenuation of the boundary perimeter within the terrace frequency band relative to the ideal state, and is the core basis for subsequent area correction.

[0046] It should be noted that the five-point difference method is used to solve for the discrete coordinates of the boundary curve. This involves taking the coordinates of two adjacent points before and after the target point and calculating the differentials of the vertical and horizontal axes of that point using the difference formula. The sensor's line-of-sight azimuth is obtained by directly extracting the sensor's orbital attitude parameters from the satellite imagery metadata file, transforming them to obtain the line-of-sight azimuth in the ground coordinate system, and simultaneously calibrating using ground control points. The calibration error must be controlled within 0.5 degrees. The initial value of the direction sensitivity coefficient is set to 1.0, and subsequent precise calibration updates can be performed using UAV data. The discretization of the boundary curve integral can be achieved by first dividing the boundary curve into segments twice the ground sampling distance, taking the parameter value at the midpoint of each segment as the representative value of that segment, and then converting the integration operation into the sum of the products of the representative values ​​of each segment and the arc length of the segment, thus realizing the numerical calculation of the integral.

[0047] In one embodiment of the present invention, the direction perpendicular to the direction of maximum slope is set as the UAV flight path layout direction, and half of the main wavelength value of the terrace frequency band is set as the flight path spacing across the contour line normal.

[0048] The key weight field is obtained by calculating the product of the terrace weight and the three-coupled perimeter gain kernel. The local maxima of the key weight field in the flight path direction are identified. Taking each local maxima as the center, the main wavelength of the terrace frequency band is extended to one-quarter of the length on both sides of the contour line normal. The resulting set of regions is determined as the key acquisition zone.

[0049] It should be noted that the drone flight path direction is defined to adapt to the characteristics of the terraced fields, specifically perpendicular to the direction of maximum slope, ensuring complete coverage of the terraced fields. The distance between drone flight paths across contour lines represents the distance between adjacent drone flight paths in the direction of the contour line, determined by the dominant wavelength of the terraced field frequency band. The key weight field represents the product of the terraced field weight and the three-coupled perimeter gain kernel, reflecting the degree of influence of the region on the perimeter gain calculation. Local maxima of the key weight field indicate points where the value of the key weight field exhibits a local peak in the flight path direction; these points correspond to highly sensitive areas of perimeter gain deviation. The key acquisition zone represents the densely acquired drone acquisition area centered on the local maxima of the key weight field.

[0050] In one embodiment of the present invention, UAV elevation data and boundary curves are acquired within a key acquisition zone. The deviation of the boundary curve segment length relative to the average length is calculated, and the deviation is multiplied by the absolute value of the cosine of the difference between the boundary tangential azimuth and the line-of-sight azimuth to obtain the micro-scale arc length jitter variance along the line of sight. Simultaneously, the root mean square of the elevation difference between the UAV elevation data and the original digital elevation model is calculated to obtain the vertical mean square error. The vertical mean square error is multiplied by the tangent of the off-axis angle and divided by the ground sampling distance to obtain an intermediate variable. The micro-scale arc length jitter variance is divided by the intermediate variable to obtain the direction sensitivity coefficient.

[0051] The error sensitivity terms of the orthophoto and elevation models are recalculated based on the direction sensitivity coefficient. These terms are then multiplied by the length scaling and visibility to obtain the updated three-coupled perimeter gain kernel. The product of the updated three-coupled perimeter gain kernel and the terrace weights is integrated along the boundary curve. The integral result is divided by the integral result of the terrace weights along the boundary curve, and one is subtracted from the quotient to obtain the updated terrace three-coupled perimeter gain.

[0052] Specifically, directional sensitivity coefficient The calculation formula is as follows:

[0053]

[0054] in This represents the variance of the microscale arc length jitter along the line of sight. This represents the vertical mean square error. Indicates the off-axis angle. Indicates the ground sampling distance.

[0055] It should be noted that UAV elevation data represents regional elevation information acquired and inverted by the measurement equipment onboard the UAV, possessing higher local accuracy than conventional digital elevation models. The boundary curves extracted by the UAV represent the outlines of the land parcels extracted based on high-resolution UAV data. The length of a boundary curve segment represents the actual arc length of each segment after discretizing the boundary curve extracted by the UAV. The average length of a boundary curve segment represents the arithmetic mean of the lengths of all discrete segments of the boundary curve, used to measure the baseline level of segment length. The deviation of the boundary curve segment length from the average length reflects the degree of fluctuation in the length of the boundary segment. The microscale arc length jitter variance along the line of sight represents the statistical variance of the boundary segment length deviation after line-of-sight azimuth correction, quantifying the intensity of jitter caused by the observation angle. The elevation difference between the UAV elevation data and the original digital elevation model represents the difference between the UAV elevation data and the original digital elevation model elevation at the same point, and is a core indicator for evaluating the error of the original elevation data. The vertical mean square error reflects the local vertical accuracy of the original digital elevation model. The direction sensitivity coefficient is obtained by the ratio of the microscale arc length jitter variance to the intermediate variable, and can calibrate the sensitivity of elevation error to plane boundary jitter. The updated orthophoto and elevation model error sensitivity terms represent the error sensitivity terms recalculated after substituting the direction sensitivity coefficient, which can improve the accuracy of observation error correction. The updated three-coupled perimeter gain kernel represents the comprehensive gain kernel that integrates the updated error sensitivity terms, which can more accurately reflect the perimeter changes caused by the superposition of multiple factors. The updated terraced three-coupled perimeter gain represents the perimeter gain calculated based on the updated gain kernel, which serves as a high-precision benchmark for subsequent area correction.

[0056] It should be noted that the updated formula for calculating the error sensitivity term of the orthophoto and elevation model is: take one, add the direction sensitivity coefficient multiplied by the original digital elevation model's vertical mean square error multiplied by the off-axis angle tangent multiplied by the product of the absolute value of the cosine of the difference between the boundary tangential azimuth and the line-of-sight azimuth, and then divide by the quotient of the ground sampling distance. Traditional UAV sampling often uses uniform grid flight paths. This invention first determines the basic flight path spacing based on the dominant wavelength of the terraced field frequency band to ensure that the terraced field strip information is not lost. Then, it uses a key weight field to lock the perimeter gain high-sensitivity area for densification, which achieves both global coverage and accurate acquisition of key areas.

[0057] It should be noted that the value of the key weight field must be greater than 1.5 times the average weight of the entire field, and the weight values ​​of the five adjacent pixels in the flight path direction must all be lower than that point to be considered a local maximum. The standard for the discrete segment length of the UAV boundary curve can be three times the ground sampling distance to segment the boundary curve, ensuring that the segment length conforms to the image resolution characteristics while effectively capturing micro-scale arc length jitter. The sampling point density for UAV elevation data within the key acquisition zone should reach one sampling point per square meter to ensure the statistical reliability of the elevation difference calculation. Furthermore, segmented data with deviations exceeding three standard deviations in the arc length jitter variance calculation are removed, and the variance is recalculated to avoid interference from outliers in the calibration results; this will not be elaborated upon here.

[0058] In one embodiment of the present invention, the perimeter of the terraced frequency band is obtained by integral of the terraced weight along the boundary curve with arc length.

[0059] Multiply the perimeter gain of the three-coupled terraces by the perimeter of the terrace frequency band, and divide the product by twice the value of pi to obtain the equivalent outer radius.

[0060] The gradient vector of the brightness field of the satellite orthophoto image on the plane is calculated. The normal brightness gradient is obtained by performing a dot product operation on the gradient vector and the normal unit vector across the contour line. At the same time, a solar direction unit vector is constructed based on the solar azimuth angle. The dot product of the normal unit vector across the contour line and the solar direction unit vector is calculated to obtain the cosine projection term. The normal brightness gradient is multiplied by the cosine projection term to obtain the weight term. The product of the weight term and the terrace weight is integrated along the boundary curve. The integral result is divided by the perimeter of the terrace frequency band to obtain the weighted average. The positive or negative sign of the weighted average is extracted as the offset sign.

[0061] It should be noted that the terrace frequency band perimeter represents the weighted perimeter obtained by integrating the terrace weights along the boundary curve of the land parcel, reflecting only the effective boundary length within the terrace frequency band. The equivalent outward expansion radius represents the equivalent geometric outward expansion scale transformed from the three-coupled perimeter gain of the terrace, enabling the mapping from perimeter changes to area changes. The satellite orthorectified image brightness field represents the two-dimensional field composed of the brightness values ​​of each pixel in the satellite orthorectified image, reflecting the illumination reflection characteristics of the land surface. The brightness field plane gradient vector represents the two-dimensional gradient of the brightness field of the satellite orthorectified image in a plane coordinate system, reflecting the rate of change of brightness in the horizontal and vertical axes, and can reflect the edge features in the image. The normal brightness gradient represents the dot product of the brightness field plane gradient vector and the unit normal vector across the contour line, reflecting the intensity of the brightness change along the contour line normal. The solar direction unit vector represents a unit length vector on the horizontal plane constructed based on the solar azimuth angle, pointing to the horizontal projection direction of the sun, used to determine the boundary offset direction in conjunction with the brightness gradient. The cosine projection term represents the dot product of the unit vector of the contour line normal and the unit vector of the sun direction, quantifying the angular relationship between the contour line normal and the horizontal direction of the sun, and is used to filter effective brightness gradient information. The weight term represents the product of the normal brightness gradient and the cosine projection term, which can integrate the correlation information between the brightness normal change and the sun's azimuth. The weighted average is the integral value of the product of the weight term and the terrace weight along the boundary curve, divided by the perimeter of the terrace frequency band; its sign reflects the overall boundary offset trend. The offset sign indicates the sign of the weighted average, taking a value of positive or negative one, used to clarify the direction of area correction, i.e., whether the boundary is offset outward or inward.

[0062] Specifically, the formula for calculating the planar gradient vector of the brightness field is to construct a two-dimensional vector by taking the first-order partial derivatives of the brightness field in the horizontal and vertical directions respectively. The formula for calculating the unit vector in the solar direction is to take the cosine value of the solar azimuth angle as the horizontal component and the sine value as the vertical component and normalize it to a unit vector. This invention establishes a correlation between the perimeter of the planar figure and the outward expansion area, and associates the perimeter gain within the terraced field frequency band with the equivalent outward expansion radius. This transforms the macroscopic relative change in perimeter into a microscopic geometric offset scale, thereby achieving precise quantification of area correction. The planar gradient of the brightness field can be calculated using the Sobel operator, that is, by performing convolution operations on the brightness values ​​in the horizontal and vertical directions using a 3×3 convolution kernel respectively, obtaining the partial derivatives in the two directions and constructing the gradient vector. The numerical discretization method for the boundary curve integral is to divide the boundary curve into segments of equal arc length at twice the ground sampling distance, take the parameter value at the midpoint of each segment as the representative value of that segment, and transform the integral operation into the sum of the product of the representative value of each segment and the arc length of the segment.

[0063] In one embodiment of the present invention, the area correction amount is obtained by multiplying the offset sign, the equivalent outer radius, and the perimeter of the terrace frequency band.

[0064] Obtain the original area calculated on the projection plane, and subtract the area correction amount from the original area to obtain the corrected area.

[0065] It should be noted that the area correction amount represents an area adjustment value determined jointly by the offset sign, the equivalent outer radius, and the perimeter of the terrace frequency band. This adjustment is used to correct area measurement deviations caused by the terrace frequency band. The original area represents the plot area directly measured on the map projection plane, without considering the perimeter gain effect related to the terrace frequency band, and serves as the baseline value for area correction. The corrected area represents the final plot area obtained by subtracting the area correction amount from the original area. This eliminates the area deviation caused by the terrace frequency band and can be directly used for land spatial planning statistics.

[0066] In one embodiment of the present invention, the square of the standard uncertainty of the equivalent outward expansion radius is calculated, and its value is equal to the sum of two terms: the first term is the square of the ratio of the terrace frequency band perimeter to twice the value of pi multiplied by the square of the standard uncertainty of the terrace three-coupling perimeter gain, and the second term is the square of the ratio of the terrace three-coupling perimeter gain to twice the value of pi multiplied by the square of the standard uncertainty of the terrace frequency band perimeter.

[0067] The square of the uncertainty of the area correction is equal to the sum of two terms: the first term is the square of the product of the offset sign and the perimeter of the terrace frequency band multiplied by the square of the standard uncertainty of the equivalent expansion radius; the second term is the square of the product of the offset sign and the equivalent expansion radius multiplied by the square of the standard uncertainty of the terrace frequency band. The square of the uncertainty of the area correction is defined as the area uncertainty.

[0068] Specifically, the square of the standard uncertainty of the equivalent outward radius. The calculation formula is as follows:

[0069]

[0070] in Indicates the perimeter of the frequency band of the terraced fields. The square of the standard uncertainty of the perimeter gain of the three-coupling terraces is expressed as... This indicates the perimeter gain of the terraced fields under three-coupling. This represents the square of the standard uncertainty of the perimeter of the terraced field frequency band.

[0071] Specifically, the square of the uncertainty of the area correction. The calculation formula is as follows:

[0072]

[0073] in Indicates the offset symbol. Indicates the perimeter of the frequency band of the terraced fields. This represents the square of the standard uncertainty of the equivalent outward expansion radius. Indicates the equivalent outward radius. This represents the square of the standard uncertainty of the perimeter of the terraced field frequency band.

[0074] It should be noted that the standard uncertainty of the three-coupling perimeter gain of terraced fields can be assessed using a Type A evaluation method. This involves performing at least ten repeated measurements on the three-coupling perimeter gain of terraced fields in the same area, calculating the experimental standard deviation of the ten sets of data, and then squared the result to obtain the standard uncertainty of the parameter. The standard uncertainty of the terraced field frequency band perimeter can be obtained by first calculating the standard uncertainty of the original boundary arc length, then, considering the normalization characteristics of the terrace weights, weighting the arc length uncertainty according to the terrace weights, and finally obtaining the standard uncertainty of the terraced field frequency band perimeter, the square of which is the corresponding statistic. The area uncertainty represents the square of the area correction uncertainty and can be used to determine whether the area data meets the accuracy requirements of the national land spatial planning. Furthermore, the ratio of the area uncertainty to the corrected area must be less than three per thousand to determine that the area data meets the accuracy acceptance requirements of the national land spatial planning.

[0075] In one embodiment of the present invention, all land parcels are divided into corresponding land category sets according to the land category identifier to which the land parcels belong, such that each land category set contains the corrected area value and area uncertainty value of all land parcels under that category;

[0076] For each land category set, the corrected area values ​​of all land parcels within the set are summed to obtain the corrected total area of ​​that land category.

[0077] For each set of land categories, calculate the square of the area uncertainty of all land parcels in the set, sum all the square values, and take the square root of the sum to obtain the total uncertainty of the land category.

[0078] Construct a final output set containing all land types, where each output item consists of a land type identifier, the corresponding corrected total area, and the corresponding total uncertainty.

[0079] It should be noted that land category identifiers are used to distinguish different land use types and reflect the planning attributes of land parcels. The total area after land category correction is the sum of the corrected areas of all land parcels within the same land category set. The total uncertainty of land categories is the square root of the sum of the squares of the area uncertainties of land parcels within the land category set. The final output set contains three core pieces of information: land category identifier, total area, and total uncertainty. The area correction processes for each land parcel within a land category are independent of each other, and their uncertainties are not significantly correlated. The method of combining the squares and square roots can accurately reflect the overall accuracy of the total area of ​​land categories, which not only meets the criteria for combining measurement uncertainties but also meets the traceability requirements for data accuracy in territorial spatial planning. Land category identifiers must adopt the first-level or second-level category codes in the current land use classification standards of territorial spatial planning, and their source must be the land parcel ownership survey data or the third national land survey / third national land survey annual update data to ensure the authority and consistency of land category classification. The final output set must be in JSON or SHP vector format. The JSON format must be organized as key-value pairs of "land category identifier-total area-total uncertainty". The SHP format must write the three pieces of information into the attribute table. The delivery carrier is an encrypted local file or an online interface of a compliant land data platform, which will not be elaborated here.

[0080] It should be noted that the interval and threshold sizes are set for ease of comparison. The size of the threshold depends on the amount of sample data and the base number set by those skilled in the art for each set of sample data, as long as it does not affect the proportional relationship between the parameter and the quantized value. Furthermore, the above formulas are all dimensionless calculations, and the formulas are derived from software simulations using a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0081] The embodiments of this example have been described above. However, this example is not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms based on the guidance of this example, and all of them are within the protection scope of this example.

Claims

1. A land area information collection system based on national land spatial planning, characterized in that, include: The data acquisition module is used to acquire digital elevation models, calculate slope, slope direction, contour line normal, generate terrace weights and terrace frequency band wavelengths through bandpass filtering, and calculate terrain elevation angle and visibility by combining solar geometric parameters. The terraced field three-coupling perimeter gain calculation module is used to segment land types in satellite orthophotos to obtain boundary curves, generate length scaling by combining slope, slope direction and projection scale factor, integrate visibility and orthophoto and elevation model error sensitive terms, construct a three-coupling perimeter gain kernel, and obtain the terraced field three-coupling perimeter gain by weighted accumulation of terraced field weights. The directional sensitivity coefficient calculation module is used to set the UAV flight path spacing according to the wavelength of the terrace frequency band, determine the key acquisition zone according to the product of the terrace weight and the three-coupled perimeter gain kernel, calculate the directional sensitivity coefficient through UAV data, and update the error sensitivity terms of the orthophoto and elevation models and the three-coupled perimeter gain of the terrace based on this. The equivalent outer radius calculation module is used to obtain the frequency band perimeter of the terraced fields through weighted accumulation of terraced field weights, determine the equivalent outer radius based on the ratio of the terraced field three-coupling perimeter gain to the perimeter, and calculate the brightness gradient based on satellite orthophoto and determine the sign based on the solar azimuth angle. The area correction module is used to take the product of the sign, the equivalent outer radius, and the perimeter of the terrace frequency band as the area correction amount, subtract the correction amount from the original area to obtain the corrected area, and at the same time calculate the area uncertainty. The summary output module is used to summarize the corrected area and corresponding uncertainty of each land type as the final output.

2. The land area information collection system based on territorial spatial planning according to claim 1, characterized in that, Calculate the first-order partial derivatives of the digital elevation model in the horizontal and vertical directions of the plane coordinate system. Determine the slope based on the arctangent of the square root of the sum of the squares of the partial derivatives in these two directions. Determine the direction of maximum slope based on the arctangent of the ratio of the partial derivatives in these two directions. Then, convert the angle value of the direction of maximum slope into a unit vector of the normal direction across the contour lines in the two-dimensional plane. Centered on each pixel, a local elevation profile is extracted along the unit vector normal to the contour line. The local elevation profile is then transformed in the frequency domain to obtain the local energy spectrum. The frequency corresponding to the maximum energy in the local energy spectrum is identified as the main peak frequency. The reciprocal of the main peak frequency is calculated to obtain the main wavelength of the terrace band. At the same time, the energy spectrum within a preset neighborhood of the main peak frequency is integrated to obtain the bandpass signal energy. The bandpass signal energy is divided by the sum of the integrals of the bandpass signal energy of the entire map to obtain the normalized terrace weight. A unit vector for the sun's direction is constructed based on the trigonometric function value of the sun's azimuth angle. The maximum elevation angle of each pixel relative to the surrounding terrain is searched and calculated along this vector to determine the terrain elevation angle. The difference between the sun's elevation angle and the terrain elevation angle is calculated, and this difference is mapped to a value between zero and one using a non-linear smoothing mapping function to obtain the visibility.

3. The land area information collection system based on territorial spatial planning according to claim 1, characterized in that, The boundary curves of the land parcels are obtained by segmenting and extracting the satellite orthophotos. The differentials of each point on the boundary curves in the vertical and horizontal directions are calculated respectively, and the arctangent value of the ratio of the two in the four quadrants is taken as the boundary tangential orientation. Calculate the square of the cosine of the difference between the boundary tangential azimuth and the maximum slope direction, multiply it by the square of the slope tangent, add one, and take the reciprocal of the square root to obtain the slope projection scaling component. Multiply the slope projection scaling component by the map projection scale factor to obtain the length scaling amount.

4. A land area information collection system based on territorial spatial planning as described in claim 3, characterized in that, Calculate the absolute value of the cosine of the difference between the boundary tangential azimuth and the sensor line-of-sight azimuth, multiply it by the direction sensitivity coefficient, the vertical mean square error of the digital elevation model, and the tangent of the off-axis angle, divide the result by the ground sampling distance, and add one to the quotient to obtain the orthophoto and elevation model error sensitivity terms. Multiply the length scaling, visibility, and orthorectification and elevation model error sensitivity terms to obtain the three-coupled perimeter gain kernel; Integrate the product of the three-coupled perimeter gain kernel and the terrace weight along the boundary curve, and simultaneously integrate the terrace weight along the boundary curve. Divide the result of the former by the result of the latter, and subtract one from the quotient to obtain the three-coupled perimeter gain of the terrace.

5. A land area information collection system based on territorial spatial planning according to claim 1, characterized in that, The direction perpendicular to the maximum slope direction is set as the drone flight path layout direction, and half of the main wavelength value of the terrace frequency band is set as the flight path spacing across the contour line normal. The key weight field is obtained by calculating the product of the terrace weight and the three-coupled perimeter gain kernel. The local maxima of the key weight field in the flight path direction are identified. Taking each local maxima as the center, the main wavelength of the terrace frequency band is extended to one-quarter of the length on both sides of the contour line normal. The resulting set of regions is determined as the key acquisition zone.

6. A land area information collection system based on territorial spatial planning according to claim 5, characterized in that, Within the key acquisition zone, UAV elevation data and boundary curves are obtained. The deviation of the boundary curve segment length relative to the average length is calculated, and this deviation is multiplied by the absolute value of the cosine of the difference between the boundary tangential azimuth and the line-of-sight azimuth. The micro-scale arc length jitter variance along the line of sight is statistically obtained. At the same time, the root mean square of the elevation difference between the UAV elevation data and the original digital elevation model is calculated to obtain the vertical mean square error. The vertical mean square error is multiplied by the tangent of the off-axis angle and divided by the ground sampling distance to obtain an intermediate variable. The micro-scale arc length jitter variance is divided by this intermediate variable to obtain the direction sensitivity coefficient. The error sensitivity terms of the orthophoto and elevation models are recalculated based on the direction sensitivity coefficient. These terms are then multiplied by the length scaling and visibility to obtain the updated three-coupled perimeter gain kernel. The product of the updated three-coupled perimeter gain kernel and the terrace weights is integrated along the boundary curve. The integral result is divided by the integral result of the terrace weights along the boundary curve, and one is subtracted from the quotient to obtain the updated terrace three-coupled perimeter gain.

7. A land area information collection system based on territorial spatial planning according to claim 1, characterized in that, The perimeter of the terrace frequency band is obtained by integrating the terrace weights along the boundary curve with arc length; Multiply the perimeter gain of the three-coupled terraces by the perimeter of the terrace frequency band, and divide the product by twice the value of pi to obtain the equivalent outer radius. Calculate the gradient vector of the brightness field of the satellite orthophoto on the plane, and perform a dot product operation between the gradient vector and the unit normal vector across the contour lines to obtain the normal brightness gradient; at the same time, construct a unit vector of the solar direction based on the solar azimuth angle, and calculate the dot product between the unit normal vector across the contour lines and the unit vector of the solar direction to obtain the cosine projection term. The weight term is obtained by multiplying the normal brightness gradient with the cosine projection term. The product of this weight term and the terrace weight is integrated along the boundary curve. The integral result is divided by the perimeter of the terrace frequency band to obtain the weighted average. The positive or negative sign of this weighted average is extracted as the offset sign.

8. A land area information collection system based on territorial spatial planning as described in claim 1, characterized in that, Multiply the offset sign, the equivalent outer radius, and the perimeter of the terrace frequency band together to obtain the area correction amount; Obtain the original area calculated on the projection plane, and subtract the area correction amount from the original area to obtain the corrected area.

9. A land area information collection system based on territorial spatial planning as described in claim 8, characterized in that, The square of the standard uncertainty of the equivalent outward expansion radius is equal to the sum of two terms: the first term is the square of the ratio of the terrace frequency band perimeter to twice the value of pi multiplied by the square of the standard uncertainty of the terrace three-coupled perimeter gain; the second term is the square of the ratio of the terrace three-coupled perimeter gain to twice the value of pi multiplied by the square of the standard uncertainty of the terrace frequency band perimeter. The square of the uncertainty of the area correction is equal to the sum of two terms: the first term is the square of the product of the offset sign and the perimeter of the terrace frequency band multiplied by the square of the standard uncertainty of the equivalent expansion radius; the second term is the square of the product of the offset sign and the equivalent expansion radius multiplied by the square of the standard uncertainty of the terrace frequency band. The square of the uncertainty of the area correction is defined as the area uncertainty.

10. A land area information collection system based on territorial spatial planning according to claim 1, characterized in that, Based on the land category identifier of the land parcel, all land parcels are divided into the corresponding land category set, so that each land category set contains the corrected area value and area uncertainty value of all land parcels under that category; For each land category set, the corrected area values ​​of all land parcels within the set are summed to obtain the corrected total area of ​​that land category. For each set of land categories, calculate the square of the area uncertainty of all land parcels in the set, sum all the square values, and take the square root of the sum to obtain the total uncertainty of the land category. Construct a final output set containing all land types, where each output item consists of a land type identifier, the corresponding corrected total area, and the corresponding total uncertainty.

Images