Road area intelligent extraction system for aerial image of unmanned aerial vehicle

Abstract

The invention relates to the technical field of image processing, in particular to a road area intelligent extraction system for unmanned aerial vehicle aerial images, and the system comprises an image noise filtering module which analyzes and inputs frequency domain information of the unmanned aerial vehicle aerial images, recognizes a target noise mode, and filters the unmanned aerial vehicle aerial images. According to the method, frequency domain noise analysis and notch filtering are combined, the sharpness of the road edge is reserved, meanwhile, high-frequency noise and periodic interference are restrained, and the visibility of a low-contrast road area is improved; local pixel spectrum dispersion calculation is introduced and is dynamically fused with morphological parameters, and direction consistency and length-width ratio constraints are combined, so that the connectivity perception capability of the fragmented non-paved road is enhanced, and the fracture misjudgment under the form-changeable terrain is reduced; a normalized vegetation index and multi-band reflectivity are combined for modeling, water reflection interference is eliminated through a blue light band standard deviation, and a space mask of a non-road surface feature interference area is constructed.

Description

Technical Field

[0001] The present invention relates to the technical field of image processing, and particularly to an intelligent extraction system for road areas in UAV aerial images. Background Art

[0002] The intelligent extraction of road areas in UAV aerial images is a solution based on image processing and computer vision, used to automatically identify and extract continuous areas of unpaved roads from high-resolution remote sensing images taken by low-altitude UAVs.

[0003] In the prior art, although morphological erosion and dilation operations can suppress noise, the fixed structure kernel size will over-smooth small roads, causing topological distortion in hilly areas. Secondly, the existing methods have a single modeling dimension for non-road interference objects, usually only using single features such as spectra or textures, lacking a multi-source feature collaborative verification mechanism. A single vegetation index is prone to confusion; and in scenarios of forest cover or temporary dirt roads, it is difficult for broken roads to achieve cross-occlusion connection based on local features, and a large amount of post-processing interpolation is required, resulting in limited efficiency and accuracy. Therefore, improvements are needed. Summary of the Invention

[0004] The purpose of the present invention is to solve the drawbacks existing in the prior art, and to propose an intelligent extraction system for road areas in UAV aerial images.

[0005] To achieve the above purpose, the present invention adopts the following technical solutions: An intelligent extraction system for road areas in UAV aerial images includes: An image noise filtering module, which analyzes the frequency domain information of the input UAV aerial image, identifies the target noise pattern, filters the UAV aerial image, and obtains the filtered image data; A spectral dispersion calculation module, based on the filtered image data, divides the UAV aerial image into multiple local region blocks, calculates the local pixel spectral dispersion value of the pixel color values within each region block; compares and judges the local pixel spectral dispersion values of each region block with a preset dispersion threshold, marks the regions lower than the preset dispersion threshold, and establishes a low spectral dispersion region map; A road connectivity analysis module, based on the low spectral dispersion region map, identifies connected components, obtains the morphological parameters of the connected regions, and merges the low spectral dispersion components according to the morphological parameters of the connected regions to generate candidate unpaved road segments; A road area constraint refinement module, based on the filtered image data, identifies vegetation-covered areas and water areas, establishes a non-road object interference area, compares the spatial positions of the candidate unpaved road segments with the non-road object interference area, removes the road segments overlapping with the non-road object interference area, and connects the remaining road segments to obtain a road network.

[0006] Preferably, the step of acquiring the filtered image data is: The pixel matrix of the input drone aerial image is converted into a frequency domain complex matrix through a two-dimensional fast Fourier transform, and the real component, imaginary component and amplitude spectrum component of the frequency domain complex matrix are extracted to generate frequency domain information including amplitude spectrum kurtosis, phase spectrum variance and frequency energy distribution density; Based on the frequency domain information, traverse the abnormal peak area in the amplitude spectrum component that exceeds 3 times the standard deviation of the average amplitude value, analyze the coordinate position, half-maximum full width and amplitude attenuation gradient of the abnormal peak area, filter the spectrum feature cluster corresponding to the noise according to the coordinate position, and identify the target noise pattern; According to the target noise pattern, a notch filter group is set at the coordinates of the abnormal peak of the amplitude spectrum component, and the center frequency of each notch filter is configured to be the horizontal and vertical coordinate values ​​of the abnormal peak, the bandwidth parameter is 1.5 times the half-maximum full width value, and the attenuation intensity is the inverse of the amplitude attenuation gradient. After filtering the frequency domain complex matrix, an inverse Fourier transform is performed to obtain filtered image data.

[0007] Preferably, the step of obtaining the local pixel spectrum dispersion value is: The filtered image data is divided into grids with 16×16 pixels as the basic unit, and the sliding step of adjacent grids is set to 8 pixels to generate a set of local area blocks; Based on the local area block set, extract the values ​​of the three channels R, G and B of all pixels in each area block in the RGB color space, calculate the mean, standard deviation, maximum and minimum values ​​of each channel, and generate a color mean vector, a color standard deviation vector, a color maximum vector and a color minimum vector of each area block; The local pixel spectrum dispersion value is calculated based on the color standard deviation vector, the color maximum value vector and the color minimum value vector.

[0008] Preferably, the step of acquiring the low-spectrum discrete region map is: Traversing the local pixel spectrum discreteness value sets of all local area blocks, comparing the local pixel spectrum discreteness value of each area block with the preset discreteness threshold one by one, determining the area blocks whose local pixel spectrum discreteness value is less than the preset discreteness threshold, and generating a preliminary marked area index set; Based on the preliminary marked area index set, spatial connectivity analysis is performed on adjacent marked areas, marked areas with a boundary distance less than 2 pixels are merged, and isolated marked areas with an area less than 10 pixels are eliminated to generate an optimized marked area index set; Based on the optimized marked area index set, the coordinates of the marked area are mapped to the original image pixel matrix, the marked area and the unmarked area are filled, and a low-spectrum discrete area map is generated.

[0009] Preferably, the steps for obtaining the morphological parameters of the connected regions are as follows: Based on the low-spectral discrete region map, traverse all the pixels in the image using the eight-neighborhood traversal algorithm, aggregate the spatially continuous labeled regions into independent connected components, extract the pixel coordinate sets and centroid coordinates of each connected component, and generate a set of connected components; Based on the set of connected components, calculate the Euclidean distances from all the pixel coordinates within each connected component to the centroid coordinate, filter out the pixels with distances less than twice the average Euclidean distance of the component to form an effective direction analysis region, and generate a set of effective direction analysis regions; Based on the set of effective direction analysis regions, calculate the direction consistency of each connected component.

[0010] Preferably, the steps for obtaining the candidate non-paved road segments are as follows: Traverse the morphological parameters of the connected regions, filter out the low-spectral discrete components according to the direction consistency, and at the same time filter out the low-spectral discrete components with an aspect ratio greater than 3:1, and determine the connected components that conform to the linear distribution trend to generate a set of candidate linear components; Based on the set of candidate linear components, calculate the distances between the endpoint coordinates of spatially adjacent components. If the distance between adjacent endpoints is less than 2 pixels and the direction angle is less than 15 degrees, merge the components into a continuous region to generate an optimized merged component set; Based on the optimized merged component set, extract the geometric center connection line of the merged region, expand the boundary of the coverage region along the connection line direction to generate a convex polygon with a uniform width, and mark it as a candidate non-paved road segment.

[0011] Preferably, the steps for obtaining the non-road ground object interference area are as follows: Based on the filtered image data, extract the reflectance values of each pixel in the near-infrared band and the red band, calculate the ratio of the near-infrared band reflectance to the red band reflectance as the preliminary vegetation index, and generate a preliminary vegetation index raster map; Based on the preliminary vegetation index raster map, statistically analyze the distribution characteristics of the vegetation index in the whole map. Set the pixels with a vegetation index greater than 0.6 and a near-infrared reflectance less than 0.3 as candidate points for the vegetation-covered area, and the pixels with a reflectance lower than 0.1 and a standard deviation of the blue band reflectance less than 0.05 as candidate points for the water body area, and generate a set of candidate points for the vegetation-covered area and a set of candidate points for the water body area; Based on the set of candidate points for the vegetation-covered area and the set of candidate points for the water body area, use morphological closing operation to merge adjacent candidate point regions, remove the isolated regions with an area less than 50 pixels, and fuse the boundaries of the vegetation and water body coverage areas to form continuous patches to generate a non-road ground object interference area.

[0012] Preferably, the steps for obtaining the road network are as follows: Based on the candidate unpaved road segments and non-road feature interference areas, extract the polygon boundary vertex coordinates of both, and use Delaunay triangulation to generate the spatial topological relationship between the road segments and the interference areas, generating a spatial topological index set; Based on the spatial topological index set, calculate the morphological difference interference score between the candidate road segments and the non-road interference areas; According to the morphological difference interference score, screen and remove the road segments, and use 5-pixel buffer analysis to connect the broken endpoints of the remaining road segments to generate a road network.

[0013] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In the present invention, by combining frequency-domain noise analysis and notch filtering, high-frequency noise and periodic interference are suppressed while retaining the sharpness of road edges, improving the visibility of low-contrast road areas; introducing the calculation of local pixel spectral dispersion and dynamically fusing it with morphological parameters, combining direction consistency and aspect ratio constraints, enhancing the connectivity perception ability of fragmented unpaved roads, and reducing false break judgments under morphologically variable terrains; using a joint model of the normalized difference vegetation index and multi-band reflectance, excluding water body reflection interference through the standard deviation of the blue light band, and constructing a spatial mask for non-road feature interference areas; based on the morphological difference scoring mechanism of the Hausdorff distance and shape similarity index, quantifying the multi-dimensional spatial relationship between road segments and interference areas, retaining candidate road segments with independent morphology and slight contact when removing overlapping areas, and reducing the over-trimming problem caused by the single index of the traditional area overlap rate. Description of the Drawings

[0014] Figure 1 It is a system flowchart of the present invention. Detailed Embodiment

[0015] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0016] Please refer to Figure 1 , the present invention provides a technical solution: An intelligent road area extraction system for UAV aerial images includes: An image noise filtering module, which analyzes the frequency-domain information of the input UAV aerial image, identifies the target noise pattern, filters the UAV aerial image, and obtains the filtered image data; The spectral dispersion calculation module divides the UAV aerial image into multiple local region blocks based on the filtered image data, calculates the local pixel spectral dispersion value of the pixel color values within each region block, compares and judges the local pixel spectral dispersion values of each region block with a preset dispersion threshold, marks the regions below the preset dispersion threshold, and establishes a low spectral dispersion region map. The road connectivity analysis module identifies connected components based on the low spectral dispersion region map, obtains the morphological parameters of the connected regions, merges the low spectral dispersion components according to the morphological parameters of the connected regions, and generates candidate unpaved road segments. The road area constraint refinement module identifies vegetation-covered areas and water areas based on the filtered image data, establishes a non-road ground object interference area, compares the spatial positions of the candidate unpaved road segments with the non-road ground object interference area, removes the road segments overlapping with the non-road ground object interference area, and connects the remaining road segments to obtain a road network.

[0017] The steps for obtaining the filtered image data are as follows: Convert the pixel matrix of the input UAV aerial image into a frequency-domain complex matrix through a two-dimensional fast Fourier transform, extract the real part component, imaginary part component, and amplitude spectrum component of the frequency-domain complex matrix, and generate frequency-domain information including amplitude spectrum kurtosis, phase spectrum variance, and frequency energy distribution density. Based on the frequency-domain information, traverse the abnormal peak regions in the amplitude spectrum component that exceed 3 times the standard deviation of the average amplitude value, analyze the coordinate positions, full width at half maximum, and amplitude attenuation gradient of the abnormal peak regions, screen the spectral feature clusters corresponding to the noise according to the coordinate positions, and identify the target noise pattern. According to the target noise pattern, set a notch filter bank at the abnormal peak coordinates of the amplitude spectrum component, configure the center frequency of each notch filter as the horizontal and vertical coordinate values of the abnormal peak, the bandwidth parameter as 1.5 times the full width at half maximum value, and the attenuation intensity as the reciprocal of the amplitude attenuation gradient. Filter the frequency-domain complex matrix and then perform an inverse Fourier transform to obtain the filtered image data.

[0018] Specifically, based on the pixel matrix of the input UAV aerial image, first call the two-dimensional fast Fourier transform algorithm to convert the pixel value matrix in the spatial domain to the frequency domain, generate a frequency-domain complex matrix of the same size but containing complex values. Then, extract the real part and the imaginary part from each complex element of this frequency-domain complex matrix to form two independent real matrices, namely the real part component matrix and the imaginary part component matrix. At the same time, according to the real part and the imaginary part corresponding to each frequency point, calculate its amplitude value. The calculation method is , thus generating an amplitude spectrum component matrix. Subsequently, based on the obtained amplitude spectrum component matrix, calculate the kurtosis of its overall pixel values to quantify the sharpness of the amplitude value distribution, and calculate the variance of the phase spectrum. The phase spectrum is obtained by calculating the arctangent function for each frequency point of the complex number Calculate the arctangent function and its variance reflects the consistency or dispersion of the phase information. Finally, based on the amplitude spectrum component matrix, calculate the energy distribution density in different frequency regions. For example, obtain the power spectrum by squaring the amplitude spectrum, and then count the total energy or average value of the low-frequency, medium-frequency, and high-frequency regions (which can be divided according to the frequency radius. For example, divide the origin radius of the frequency domain into , , three regions, where and are radius thresholds preset according to the image size and noise characteristics. For example, for a 1024x1024 image, can be set to 50, can be set to 200), and summarize the three indicators of the kurtosis of the amplitude spectrum, the variance of the phase spectrum, and the frequency energy distribution density to obtain the frequency domain information characterizing the image noise characteristics.

[0019] Based on the frequency domain information generated in the previous step, especially the amplitude spectrum component matrix, first calculate the arithmetic mean and standard deviation of all amplitude values in this matrix, and set an abnormal peak determination threshold. The calculation method of this threshold is the average amplitude value plus 3 times the standard deviation. For example, if the calculated average amplitude value of the entire amplitude spectrum component is 55.2 and the standard deviation is 8.5, then the abnormal peak determination threshold is set to . Then, traverse each pixel point (frequency point) in the amplitude spectrum component matrix, compare its amplitude value with the calculated abnormal peak determination threshold of 80.7. If the amplitude value of a certain point is greater than 80.7, mark it as an abnormal point. Apply the eight-neighborhood connectivity algorithm to all marked abnormal points to aggregate the spatially adjacent abnormal points to form independent abnormal peak regions. Then, for each identified abnormal peak region, analyze its key attributes, including determining the coordinates of the point with the largest amplitude value in this region as the peak center coordinates, calculate the full width at half maximum (FWHM) of this peak along the horizontal and vertical directions at the peak center , especially their distribution patterns in the frequency domain, such as whether there are pairs of peak points symmetric about the origin, whether they are arranged on a specific straight line and other periodic noise characteristics, to screen and summarize the spectral feature clusters representing specific interference sources, and obtain the identified target noise patterns.

[0020] According to the specific information (coordinates, full width at half maximum, amplitude attenuation gradient) of each abnormal peak contained in the target noise pattern identified in the previous step, and the frequency-domain complex matrix generated in the first step, for each abnormal peak in the target noise pattern, design a corresponding two-dimensional notch filter, and construct a set of notch filters. When specifically configuring, set the center frequency of each notch filter to the center coordinate value, and set the bandwidth parameter of the filter according to the full width at half maximum (FWHM) of the peak, specifically set to 1.5 times the measured value of the full width at half maximum. For example, if the full width at half maximum of an abnormal peak is measured as 4 frequency units in the direction and 6 frequency units in the direction, then the bandwidth parameter of the corresponding notch filter can be set to cover units in the direction and units in the direction. The attenuation intensity of the filter is determined according to the reciprocal of the amplitude attenuation gradient of the peak. For example, if the calculated amplitude attenuation gradient is 0.8, then the attenuation intensity-related parameter is set to This value will be used to control the attenuation depth of the notch filter at the center frequency. Combine all the configured individual notch filters into a total filter transfer function

[0021] The steps for obtaining the local pixel spectral dispersion value are as follows: Divide the filtered image data into grids with 16×16 pixels as the basic unit, set the sliding step of adjacent grids to 8 pixels, and generate a set of local region blocks; Based on the set of local region blocks, extract the values of the R, G, and B channels of all pixels in each region block in the RGB color space, calculate the mean, standard deviation, maximum value, and minimum value of each channel, and generate the color mean vector, color standard deviation vector, color maximum value vector, and color minimum value vector of each region block; Based on the color standard deviation vector, color maximum vector, and color minimum vector, calculate the local pixel spectral dispersion value. The calculation formula is as follows: ; where, is the local pixel spectral dispersion value of the s-th regional block, is the standard deviation of color channel c, is the range of color channel c, is the absolute value of the covariance between the current channel c and the mean of the other two channels.

[0022] Specifically, based on the filtered image data obtained in the previous step, perform sliding window block processing. Specifically, set the basic unit size to 16 pixels by 16 pixels as the size of each local regional block. Traverse the entire filtered image data. Starting from the top-left coordinate (0, 0) of the image, extract the first 16×16 pixel block. Subsequently, set the sliding step sizes in both the horizontal and vertical directions to 8 pixels. This means that the next horizontally adjacent block will be extracted starting from the coordinate (8, 0) with 16×16 pixels, and the next one will start from (16, 0), and so on, until reaching near the right boundary of the image (the starting column coordinate of the last block satisfies is less than the width of the image). Similarly, after completing the extraction of the first row of blocks, move to the next row. The starting block coordinate is (0, 8), and the subsequent blocks are (8, 8), (16, 8), etc., until covering the entire image area (the starting row coordinate of the last block satisfies is less than the height of the image). This sliding window method with overlap ensures a smooth transition of pixel neighborhood information between blocks. Organize the coordinates and corresponding pixel data of all the extracted 16×16 pixel blocks to generate a local regional block set.

[0023] Based on the local regional block set generated in the previous step, perform color feature statistics on each local regional block in the set (i.e., a 16×16 pixel block). First, access all 256 pixel points within the regional block and read the three-channel values of each pixel point in the RGB color space, namely the red (R) channel value, green (G) channel value, and blue (B) channel value. These values are usually in the integer range from 0 to 255. Then, for the 256 pixels within the regional block, perform statistical calculations on the numerical lists of the R, G, and B channels respectively. Calculate the arithmetic mean of each channel's numerical values to obtain the average color of the block, calculate the standard deviation of each channel's numerical values to quantify the dispersion degree of the color values within the channel, find and record the maximum and minimum values in each channel's numerical values. Subsequently, combine the average values of the three channels calculated into a three-dimensional vector, namely the color mean vector (for example ), combine the standard deviations of the three channels into a color standard deviation vector (e.g., ), combine the maximum values of the three channels into a color maximum value vector (e.g., ), and combine the minimum values of the three channels into a color minimum value vector (e.g., ). Repeat this process for all local region blocks, and finally generate a color mean vector, a color standard deviation vector, a color maximum value vector, and a color minimum value vector corresponding to each region block.

[0024] Formula: , the benefit of the formula is that it comprehensively considers the variations within each color channel (standard deviation and range ) and the inter-channel relationships (covariance ), thus more accurately quantifying the spectral uniformity or complexity of local regions. The standard deviation reflects the overall fluctuation of pixel values, the range reflects the dynamic range of colors, and the covariance term adjusts the influence of the range by examining the correlation between the current channel and the means of the other two channels. When the color channels within a region co-vary (e.g., an overall increase or decrease in brightness causes synchronous increases or decreases in R, G, and B, with a large absolute value of covariance), it indicates that its internal structure may be relatively simple (such as a homogeneous surface under shadow or lighting changes). At this time, even if the range is large, through the increase of the denominator , the contribution of the range to the total dispersion value will be reduced. Conversely, if the variation relationships between channels are chaotic (small absolute value of covariance), the greater influence of the range is retained, which can better reflect the characteristics of mixed pixels or complex texture regions. This design makes give lower values to regions such as roads with relatively single spectra and relatively stable internal color relationships, and higher values to complex regions such as vegetation and building edges, improving the discrimination; The steps to obtain parameter are as follows: This parameter represents the standard deviation of the pixel values of color channel within the -th region block ( is one of R, G, B), and it directly comes from the corresponding component in the color standard deviation vector calculated for each local region block. For example, for the -th region block, directly extract it from its color standard deviation vector . If the vector value is , then , , .

[0025] The steps to obtain parameter are as follows: This parameter represents the color channel within the The range of pixel values, that is, the difference between the maximum and minimum values of this channel, is the color maximum value vector obtained through calculation and the color minimum value vector are calculated. For example, for the th region block, its color maximum value vector is , and the color minimum value vector is , then the range of the R channel , the range of the G channel , and the range of the B channel .

[0026] Parameter is obtained as follows: This parameter represents the covariance between the list of pixel values of the current color channel within the th region block and the corresponding mean value list of the pixel values of the other two channels. For example, the covariance of the R channel of the th block with the mean value lists of its G and B channels is , the covariance of the G channel with the mean value lists of its R and B channels is , and the covariance of the B channel with the mean value lists of its R and G channels is , then the absolute values used in the formula are respectively , , .

[0027] Calculation process: Taking the th local region block as an example, substitute the parameter values obtained previously for calculation: Known parameters: , , , , , , ; Calculate the contribution of the R channel: ; Calculate the contribution of the G channel: ; Calculate the contribution of the B channel: ; Calculate the total local pixel spectral dispersion value : ; This result shows that the calculated result of the local pixel spectral dispersion value of the th local region block is , and this value comprehensively reflects the degree of color change within the 16×16 pixel region and the correlation between channels. A lower value (for example, less than a preset dispersion threshold, and this threshold needs to be based on the road and non-road regions in a large number of sample images The distribution statistics are determined, and the boundary value that can better distinguish the two is selected, such as 30. It tends to indicate that the spectral characteristics of the area are relatively simple and stable, and it may be part of a homogeneous area such as roads and bare land. A low value (e.g. greater than 30) indicates that the spectral complexity of the area is high, which may be a reflection of vegetation, buildings, mixed land features, etc.

[0028] The steps for obtaining the low-spectral discrete area map are: Traversing the local pixel spectrum discreteness value sets of all local area blocks, comparing the local pixel spectrum discreteness value of each area block with the preset discreteness threshold one by one, determining the area blocks whose local pixel spectrum discreteness value is less than the preset discreteness threshold, and generating a preliminary marked area index set; Based on the preliminary marked area index set, the spatial connectivity analysis of adjacent marked areas is performed, the marked areas with boundary distance less than 2 pixels are merged, and the isolated marked areas with an area less than 10 pixels are eliminated to generate an optimized marked area index set; Based on the optimized marked area index set, the coordinates of the marked area are mapped to the original image pixel matrix, and the marked area and the unmarked area are filled to generate a low-spectral discrete area map.

[0029] Specifically, the local pixel spectral discreteness values ​​of all local area blocks obtained in the previous step are traversed The set contains the discrete values ​​calculated for each 16×16 pixel area block. The value is different from a pre-set discreteness threshold For comparison, the preset discreteness threshold The determination process is as follows: select a representative set of drone aerial images covering a variety of scenes (different lighting, weather, and surface cover), manually and accurately outline typical road areas (especially unpaved roads) and non-road areas (such as vegetation, buildings, water bodies, bare land, etc.) in the sample images, and apply the same 16×16 overlapping blocks and Calculation process, get a large number of known categories Value, statistical road area Values ​​and non-road areas The distribution histogram of the values ​​​​is used to analyze the overlap of the two types of distributions and select a value that can better distinguish the two types of areas. The value is used as the threshold. For example, the maximum inter-class variance method can be used to automatically calculate a threshold, or a value between the peaks of the two distributions can be selected based on empirical observations. For example, if statistical analysis shows that most road area blocks have The values ​​are concentrated between 15 and 35, but not the road area blocks. If most of the values are distributed above 40, the preset dispersion threshold can be set to 30. During the traversal process, for the nth region block, if its value is less than (for example ), then it is determined that this region block belongs to the low-spectral dispersion region, and the index of this region block (such as its row, column number or starting pixel coordinates in the grid) is recorded. The indices of all region blocks that meet the conditions are gathered to generate a preliminary marked region index set.

[0030] Based on the preliminary marked region index set generated in the previous step, optimize the spatial relationship of these region blocks preliminarily marked as low-spectral dispersion. First, map the position information of these marked region blocks to a two-dimensional space (corresponding to the original image space), regard each marked block as a spatial entity, conduct a spatial connectivity analysis, and find marked region blocks that are adjacent or close in space. Specifically, calculate the minimum boundary distance between any two independent and marked connected regions (initially each marked block is an independent connected region). This distance is obtained by calculating the Euclidean distance from all pixel points on the boundary of one region to all pixel points on the boundary of another region and taking the minimum value. Set a merging distance threshold This threshold is set according to experience and is used to connect regions belonging to the same road that are separated by minor discontinuities (such as small cracks on the road surface, short-term occlusions or image noise). For example, set to 2 pixels. If the minimum boundary distance between two independent marked regions is less than (for example, the calculated distance is 1.5 pixels, less than 2 pixels), then these two regions are merged into a larger connected region. Repeat this merging process until there are no more regions that meet the conditions for merging. Next, perform area filtering on all connected regions formed after the merging process. Calculate the total pixel area covered by each connected region (counting the pixel points in all 16x16 blocks that make up the region after removing duplicates), and set a minimum area threshold , which is used to eliminate low-dispersion regions that are too small and are likely to be noise or accidentally formed by non-road features. This threshold is empirically determined based on the minimum effective width and length of the road and the image resolution in actual applications. For example, for an application aiming to extract the carriageway, can be set to 10 pixels. If the total pixel area of a certain connected region is less than (for example, the calculated area is 8 pixels, less than 10 pixels), then it is determined that this region is an isolated noise or invalid region and is removed from the marked regions. Retain all connected regions with an area greater than or equal to . The corresponding region block indices form an optimized marked region index set.

[0031] Based on the optimized marked region index set generated in the previous step, which contains the region block indexes considered as valid low-spectral discrete regions after merging and area filtering, perform mapping and filling operations to generate the final binary image, that is, the low-spectral discrete region map. The specific process is as follows: Create a blank binary image matrix with exactly the same size as the original filtered image data (for example pixels), and initialize all pixel values to the background value, usually set to 0, representing non-marked regions. Then, traverse each region block index in the optimized marked region index set, and determine the specific pixel coordinate range of this 16×16 region block in the original image coordinate system according to the index (for example, the block index corresponds to the top-left pixel coordinate , and the coverage range is ). Next, modify the pixel values at the corresponding positions of all pixel points (a total of 256) covered by this region block in the blank binary image matrix to the foreground value, usually set to 1, representing the marked region. Since there is an 8-pixel overlap between region blocks, some pixel points may be covered by multiple marked region blocks and be assigned the value 1 multiple times, which does not affect the final result. When all indexes in the optimized marked region index set have been traversed, the binary image matrix is constructed. The region with pixel value 1 represents the finally determined low-spectral discrete region, and the region with pixel value 0 represents other regions, thus obtaining the low-spectral discrete region map.

[0032] The steps for obtaining the morphological parameters of connected regions are as follows: Based on the low-spectral discrete region map, traverse all pixels in the image using the eight-neighborhood traversal algorithm, aggregate spatially continuous marked regions into independent connected components, extract the pixel coordinate set and centroid coordinate of each connected component, and generate a connected component set; Based on the connected component set, calculate the Euclidean distance from all pixel coordinates within each connected component to the centroid coordinate, and filter out the pixels with a distance less than twice the average Euclidean distance of the component to form an effective direction analysis region, and generate an effective direction analysis region set; Based on the effective direction analysis region set, calculate the direction consistency of each connected component. The calculation formula is: ; where, is the direction consistency of the i-th connected component, is the total number of pixels in the effective direction analysis region of the i-th component, is the azimuth angle of the j-th pixel in the i-th component relative to the centroid, and are the horizontal and vertical coordinates of the j-th pixel in the i-th component, and are the horizontal and vertical coordinates of the centroid of the i-th component.

[0033] Specifically, based on the low-spectral discrete region map generated in the previous step (where marked pixels, e.g., with a value of 1, represent low-discrete regions, and unmarked pixels, e.g., with a value of 0, represent other regions), each independent marked region in the map is identified and separated. Specifically, the eight-neighborhood traversal algorithm is adopted. This algorithm scans row by row starting from the first pixel of the image (e.g., the upper left corner). When an unvisited marked pixel (with a value of 1) is encountered, it is used as the starting point (seed point) of a new connected component, and a unique identifier is assigned to this component. At the same time, the pixel coordinate list, the sum of coordinates (used to calculate the centroid), and the pixel counter of this component are initialized. Then, a breadth-first search or depth-first search is performed using a queue or stack structure. The seed point is added to the queue / stack and marked as visited. The pixels in the queue / stack are processed in a loop: take out a pixel, check its eight neighboring pixels (up, down, left, right, upper left, upper right, lower left, lower right). For each neighboring pixel, if it is within the image range, is a marked pixel (with a value of 1), and has not been visited, it is marked as visited, given the same component identifier, its coordinates are added to the pixel coordinate list of the current component, its coordinate values are accumulated into the sum, the pixel counter is incremented, and it is added to the queue / stack. When the queue / stack is empty, it indicates that all marked pixels connected to the seed point have been found, forming a complete independent connected component. Record the set of pixel coordinates of this component (a list containing all the pixels of this component and the centroid coordinates calculated by dividing the sum of coordinates by the pixel counter , that is and , where is the total number of pixels in the component. Continue to scan the image until all marked pixels are visited and belong to a certain connected component, and finally generate a set of connected components containing information about all independent connected components.

[0034] Based on the set of connected components generated in the previous step, where each component contains its set of pixel coordinates and centroid coordinates, further screen the pixels inside each connected component to extract the region that can better represent its main shape and direction. For the th connected component, first traverse each pixel in its set of pixel coordinates , calculate the Euclidean distance from this pixel to the centroid of this component. The calculation method is . Repeat this calculation for all pixels in this component to obtain a set of distance values. Then, calculate the arithmetic mean of these distance values to obtain the average Euclidean distance of this component. Next, set a screening threshold, which is 2 times the average Euclidean distance of the component, that is , this multiple (2 times) is set based on experience, aiming to retain the pixels in the core area of the component while removing those pixels that are far from the centroid and may belong to branches or noise. These points far from the center may interfere with the subsequent judgment of the main direction of the component. For example, if the average Euclidean distance of a connected component is calculated to be 12.5 pixels, then the screening threshold is pixels. Finally, traverse all the pixels of this component again , compare the calculated Euclidean distance with the threshold . If is less than (for example ), then retain this pixel and regard it as a valid pixel. All the retained pixels together constitute the effective direction analysis area of the th connected component. The same operation is performed on all connected components to generate a set of effective direction analysis areas.

[0035] Formula: , the benefit of the formula is that it provides a method to measure the consistency of the distribution direction of pixel points inside the connected component, and is especially suitable for judging whether the extracted low-spectral discrete area has the linear or slender shape usually presented by roads. The formula calculates the azimuth angle of each pixel relative to the centroid and doubles it to , solving the 180-degree ambiguity problem of direction (that is, no matter whether the pixel is on one side of the centroid or on the completely opposite side, as long as it is on the same straight line, its value is the same in the sense of modulo ), which enables the formula to give a high consistency score for shapes that extend bidirectionally along a single axis (such as straight line segments); The steps to obtain the parameter are as follows: This parameter represents the total number of pixels included in the effective direction analysis area of the th connected component. For example, if the number of remaining valid pixels after screening the th connected component is 550, then .

[0036] The steps to obtain the parameter are as follows: This parameter represents the azimuth angle of the th pixel in the effective direction analysis area of the th connected component relative to the centroid of this component. For example, the centroid of the th component is , and the The pixel coordinates are , then the azimuth angle of this pixel radians.

[0037] Parameter and are obtained as follows: These two parameters respectively represent the th abscissa and ordinate of the th pixel in the effective direction analysis region of the connected component. These coordinates are directly from the pixel list of the effective direction analysis region, and this list contains the original image coordinates of all pixels that meet the distance screening conditions. For example, for the th pixel in the above example, its coordinates are , .

[0038] Parameter and are obtained as follows: These two parameters respectively represent the abscissa and ordinate of the centroid of the th connected component. For example, the centroid coordinates of the th component are , .

[0039] Calculation process: Taking the th connected component as an example, its effective direction analysis region contains pixels, the centroid is , the pixel coordinates and the calculated azimuth angle and related values are as follows: Pixel : , , rad, rad, , . Pixel : , , rad, rad( ), , . Pixel : , , rad, rad, , .

[0040] Calculate : ; Calculate : ; calculate : ; This result shows that for this simplified example connected component containing 3 valid pixels , its directional consistency The calculated value is , this value is between 0 and 1, The closer the value is to 1, the more the pixels that make up the effective area tend to be distributed on a straight line passing through the centroid, showing a strong linear feature. This is usually the morphological manifestation of slender features such as roads. The closer the value is to 0, the more diffuse the pixel distribution is, lacking obvious directionality, and more like a clumpy area.

[0041] The steps for obtaining candidate unpaved road segments are as follows: Traverse the morphological parameters of the connected region, select low-spectral discrete components according to directional consistency, and at the same time select low-spectral discrete components with an aspect ratio greater than 3:1, determine the connected components that meet the linear distribution trend, and generate a set of candidate linear components; Based on the candidate linear component set, the distance between the endpoint coordinates of the spatially adjacent components is calculated. If the distance between the adjacent endpoints is less than 2 pixels and the direction angle is less than 15 degrees, the merged component is a continuous area, and an optimized merged component set is generated; Based on the optimized merge component set, the geometric center line of the merged area is extracted, and the boundary of the covered area is extended along the line direction to generate a convex polygon with uniform width, which is marked as a candidate unpaved road segment.

[0042] Specifically, traverse each connected component obtained in the previous step and its associated morphological parameters, especially the directional consistency Value, filter these low spectral discrete components, set a direction consistency threshold , which is determined based on prior knowledge or statistical analysis of sample data, and is intended to retain components with obvious linear features, for example, by analyzing the known road sample and non-road sample components. Value distribution, find the road component The value is usually higher than 0.7, while non-road components (such as patchy farmland and parking lots) The value is low, so you can set If a component Value greater than (For example ), it is preliminarily considered that it meets the requirements of linear characteristics. At the same time, the aspect ratio of the component is calculated. First, the principal component analysis method is used to calculate the major axis direction and minor axis direction of the component pixel coordinate set, and then the maximum projection length of the component in the major axis direction is calculated. and the maximum projection length in the minor axis direction , the aspect ratio is . Set an aspect ratio threshold . According to the characteristic that roads are usually slender in shape, this threshold is set to 3:1, that is . If the aspect ratio of the component is greater than (for example, calculated as ), then it is considered that its shape conforms to the slender feature. Only when the direction consistency is greater than and the aspect ratio is greater than for the connected components that meet these two conditions, can they be determined to conform to the linear distribution trend and be selected and added to the candidate linear component set.

[0043] Based on the candidate linear component set generated in the previous step, further merge and connect these linear components that may belong to the same road but are separated. First, it is necessary to determine the two endpoint coordinates of each candidate linear component, which can be achieved by finding the two pixel points with the farthest projection distance along the main axis direction of the component (obtained when calculating the aspect ratio), or by extracting the two endpoints of the skeleton line after skeletonizing the component (such as using morphological thinning). At the same time, determine the main direction vector of each component, for example, represented by the vector connecting its two endpoints or the main axis direction vector. Next, traverse each pair of candidate linear components that may be adjacent in space in the set (for example, the distance between their minimum bounding rectangles is less than a preset large search radius, such as 50 pixels), calculate the Euclidean distance between the closest pair of endpoints between them , and at the same time calculate the angle between the main direction vectors of these two components (take the acute or obtuse angle between 0 and 180 degrees). Set two merge condition thresholds: the endpoint distance threshold and the direction angle threshold . These two thresholds are set according to experience to allow connecting broken road segments with small gaps and basically the same direction. For example, set pixels, degrees. If the shortest endpoint distance of a pair of adjacent components is less than (for example ) and their direction angle is less than (for example ), then it is determined that these two components should be merged, merge their pixel sets, and update the relevant attributes of the merged area (such as recalculating the endpoints and the main direction). Repeat this process, using, for example, a graph-based connection method, until all adjacent components that meet the conditions are merged, generating an optimized merged component set.

[0044] Based on the optimized merged component set generated in the previous step, where each element represents a merged and more continuous potential road area, perform geometric normalization on each merged area. First, extract the connecting line of the geometric centers of the merged area, which can be regarded as the axis of the road segment. The extraction method can be the main axis segment obtained by calculating the principal component analysis of the pixel point set of the area, or performing skeletonization on the area to extract the longest skeleton path, or simply connecting the centroids of the original (before merging) candidate linear components included in the merged area and performing smoothing to obtain a polyline segment. After obtaining the center line, set a road width parameter , which is set empirically according to the typical width of unpaved roads and the image resolution in the study area. For example, if the target road width is about 5 meters and the image resolution is 0.5 meters / pixel, then it can be set pixels. Then, along each point of the extracted connecting line of the geometric centers, extend in both directions perpendicular to it by distance (i.e., extend 5 pixels in each direction), forming a series of equal-width line segments perpendicular to the center line. All the endpoints of these line segments outline the set of boundary contour points of the road area. Finally, apply the convex hull algorithm to this set of boundary contour points to calculate the smallest convex polygon that contains all the boundary points. This convex polygon has the characteristics of approximately uniform width (controlled by ) and regular shape, and mark this generated convex polygon as a candidate unpaved road segment.

[0045] The steps for obtaining the non-road feature interference area are as follows: Based on the filtered image data, extract the reflectance values of each pixel in the near-infrared band and the red band, calculate the ratio of the near-infrared band reflectance to the red band reflectance as the preliminary vegetation index, and generate a preliminary vegetation index raster map; Based on the preliminary vegetation index raster map, statistically analyze the distribution characteristics of the vegetation index of the whole map. Set the pixels with a vegetation index greater than 0.6 and a near-infrared reflectance less than 0.3 as candidate points for the vegetation-covered area, and the pixels with a reflectance lower than 0.1 and a standard deviation of the blue band reflectance less than 0.05 as candidate points for the water area, and generate a set of candidate points for the vegetation-covered area and a set of candidate points for the water area; Based on the set of candidate points for the vegetation-covered area and the set of candidate points for the water area, use morphological closing operation to merge adjacent candidate point areas, remove isolated areas with an area less than 50 pixels, and fuse the boundaries of the vegetation and water-covered areas to form continuous patches, generating the non-road feature interference area.

[0046] Specifically, based on, for example, the available filtered image data containing the near-infrared and red bands (the filtered image data here refers to the multi-spectral or hyperspectral image data that has undergone preprocessing, such as radiometric calibration and atmospheric correction, to obtain the surface reflectance values), perform operations on each pixel point in the image. First, for the pixel position , extract their reflectance values in the near-infrared band respectively and the reflectance values in the red light band . These reflectance values are dimensionless and usually range from 0 to 1. Subsequently, calculate the ratio of these two reflectance values, that is . During the calculation process, if the red band reflectance is extremely small or zero, to avoid division by zero error, a small positive number can be set to replace it (for example, 0.0001) or the result can be set to a specific extremely large value or an invalid value mark. This ratio , as a simplified preliminary vegetation index, the magnitude of its value is usually positively correlated with the density of surface vegetation cover. Because healthy vegetation strongly reflects near-infrared light and strongly absorbs red light, store the values calculated for all pixels in the image in a new two-dimensional raster data structure, which has the same spatial resolution and range as the original image, to generate a preliminary vegetation index raster map

[0047] . Based on the preliminary vegetation index raster map generated in the previous step and combined with the original filtered image data containing reflectances in bands such as near-infrared, red, and blue, classify pixels to identify candidate areas for vegetation and water bodies. First, statistically analyze the distribution characteristics of all pixel values in the preliminary vegetation index raster map, such as calculating its histogram, mean, and standard deviation, to understand the approximate range of the vegetation index in the study area. Set the determination criteria for candidate points in the vegetation-covered area: the preliminary vegetation index needs to be greater than the threshold , and its near-infrared band reflectance needs to be less than the threshold . The threshold is set according to experience. For example, by comparing the index values of known vegetation areas, it is set to 0.6, indicating that the NIR reflectance is significantly higher than the Red reflectance; the threshold is also set according to experience. For example, it is set to 0.3 to exclude some non-vegetation ground objects with high reflectance (such as the roofs of some buildings) or specific types (such as very dense and dark-colored) of vegetation. For candidate points in the water body area, set the determination criteria: its reflectance needs to be lower than the threshold , and the standard deviation of its blue band reflectance in the local neighborhood needs to be less than the threshold . Here , the blue band reflectance can be selected because the water body usually has the lowest reflectance in this band. The threshold is set to 0.1 according to experience, reflecting the low reflectance characteristics of the water body; the local neighborhood standard deviation It is obtained by calculating the standard deviation of the reflectance of all pixels in the blue light band within a window of, for example, 5×5 around each pixel. This value is used to measure the texture uniformity of the area. The water surface is usually very uniform, so a relatively small value, such as 0.05, is set to distinguish it from shadows or asphalt roads with the same low reflectance but possible texture. Each pixel in the image is traversed and judged according to the above criteria: If the vegetation condition is met (for example, and ), then its coordinates are recorded as vegetation coverage candidate points; if the water body condition is met (for example, and ), then its coordinates are recorded as water body area candidate points, and finally a set of vegetation coverage candidate points and a set of water body candidate points are generated.

[0048] Based on the set of vegetation coverage candidate points and the set of water body candidate points generated in the previous step, both of which are lists of pixel coordinates, morphological processing and spatial aggregation are performed on these discrete candidate points. First, the set of vegetation candidate points and the set of water body candidate points are respectively mapped onto two independent binary images. The candidate point positions are set to the foreground value 1, and the rest are set to the background value 0. Morphological closing operations are independently performed on these two binary images. The closing operation consists of dilation followed by erosion. A small structuring element is selected, such as a 3×3 pixel-sized all-1 square kernel. The size of this structuring element is selected according to experience and is designed to fill small gaps between candidate points and small holes inside the area, so that originally close but disconnected candidate point areas can be connected to form more complete patches. After the closing operation, the connected component analysis (using eight-neighbor connection) is performed on the two binary images respectively to identify all independent connected regions composed of foreground pixels, and the pixel area of each connected region is calculated. A minimum area threshold is set. This threshold is set according to experience and is used to filter out those small-area, likely noise or insignificant sporadic patches. For example, 50 pixels are set. If the area of a connected region is less than 50 pixels, then all pixel values within this region are set from 1 to 0 and removed, and all regions with an area greater than or equal to 50 pixels are retained. Finally, a logical OR operation is performed on the binary image of the vegetation area and the binary image of the water body area after area filtering, that is, a new blank binary image is created. If a pixel is at least 1 (foreground value) in the processed vegetation map or the processed water body map, then this pixel is set to 1 in the new map, otherwise it is set to 0. In this way, the verified vegetation area and water body area are fused into continuous patches to generate the final non-road ground object interference area.

[0049] The steps for obtaining the road network are as follows: Based on the candidate unpaved road segments and non-road object interference areas, extract the polygon boundary vertex coordinates of both, use Delaunay triangulation to generate the spatial topological relationship between the road segments and the interference areas, and generate a spatial topological index set; Based on the spatial topological index set, calculate the morphological difference interference score between the candidate road segments and the non-road interference areas. The calculation formula is: ; Among them, is the morphological difference interference score of the k-th road segment; , is the boundary point set of the k-th road segment and the boundary point set of the m-th interference area The Hausdorff distance of; , is the shape similarity index, and are the perimeters of the road segment and the interference area respectively, is the Pixel area of the candidate unpaved road segment, is the Pixel area of the non-road object interference area; is the smoothing coefficient to prevent the denominator from being zero; is the overlapping area between the road segment and the interference area; is the total number of non-road object interference areas; According to the morphological difference interference score, screen and remove the road segments, and use a 5-pixel buffer analysis to connect the broken endpoints of the remaining road segments to generate a road network.

[0050] Specifically, based on the candidate unpaved road segments and non-road object interference areas obtained in the previous steps, it is first necessary to clarify the spatial proximity relationship between the two, extract the boundary vertex coordinate sequences of each candidate unpaved road segment polygon and each non-road object interference area polygon, and merge all these vertex coordinates from the road segments and interference areas into a unified vertex set. Perform a Delaunay triangulation operation on this set containing all the key geometric feature points (boundary vertices). This algorithm, such as the Bowyer-Watson algorithm, will generate a triangular mesh covering the entire area, and the vertices of the mesh are the input boundary vertices. This triangular mesh has the empty circumcircle property, that is, the interior of the circumcircle of any triangle does not contain any other vertices. By checking the edges of the generated triangular mesh, it is possible to effectively identify which road segment vertices are directly connected (sharing a triangular mesh edge) or very close (belonging to the same or adjacent triangles) to the interference area vertices. This adjacency relationship based on shared vertices or edges defines the spatial topological relationship between the road segments and the interference areas. Record all pairs of (road segment index, interference area index) with such topological proximity relationships to generate a spatial topological index set.

[0051] Formula: , the benefit of the formula is that it provides a quantitative index for evaluating the degree to which each candidate unpaved road segment is affected by the non-road feature interference area . This index comprehensively considers multiple factors: First, the overlapping area in the numerator directly reflects the degree of spatial conflict between the road segment and the interference area. The larger the overlap, the more serious the potential interference; Second, the in the denominator (Hausdorff distance) measures the proximity of the boundaries of the two. The farther the distance, even if there is an overlap, the influence weight will be reduced; The in the denominator (shape similarity index) compares the shape compactness (perimeter / area ratio) of the two. If the shapes of the two are very different (e.g., a slender road and a circular interference area), even if the distance is close and there is an overlap, it may indicate that the interference area has little association with the road structure, and by increasing to reduce its influence weight; ; The acquisition steps of parameter are as follows: This parameter represents the pixel area of the th candidate unpaved road segment. After the candidate unpaved road segment is generated in polygon form, it is obtained by calculating the number of pixels covered by the polygon. For example, after calculation, the th candidate road segment polygon covers 1500 pixels, then .

[0052] The acquisition steps of parameter are as follows: This parameter represents the pixel area of the th non-road feature interference area. The acquisition method is similar to that of . If the non-road feature interference area is represented by a polygon, calculate the number of pixels it covers; if it is represented by a raster image (binary image), directly count the total number of pixels with a value of 1 (representing the interference area). This value reflects the size of the interference area. For example, the area of the th interference area (such as a piece of vegetation) is calculated as 6200 pixels, then .

[0053] The acquisition steps of parameter are as follows: This parameter represents the overlapping pixel area between the th candidate unpaved road segment and the th non-road feature interference area. It is necessary to calculate the spatial intersection of the two areas (polygon or raster). If both are polygons, the intersection polygon can be obtained using polygon clipping or Boolean operations in computational geometry, and then calculate its area; For example, it is calculated that the overlapping area between the th road segment and the If 450 pixels overlap in an interference area, then .

[0054] The parameter is obtained as follows: This parameter represents the th set of boundary points of a road segment and the th set of boundary points of an interference area . First, extract and as a list of discrete point coordinates, then calculate the standard Hausdorff distance . Finally, , take the larger value of it and 1.0 to avoid the influence of a distance of 0 or too small a value on the denominator. For example, if the calculated standard Hausdorff distance between and is 8.2 pixels, then .

[0055] The parameter is obtained as follows: This parameter represents the similarity index of the shape (compactness) between the th road segment and the th interference area. First, calculate their perimeters and respectively (obtained by accumulating the lengths of the line segments between the boundary vertices) and their areas and as before. Then, calculate their perimeter - area ratios and . Finally, take the absolute value of the difference between these two ratios, that is, . The smaller this value is, the more similar the compactness of their shapes. For example, the perimeter of the th road segment is pixels and its area is pixels; the perimeter of the th interference area is pixels and its area is pixels, then .

[0056] The parameter is obtained as follows: This parameter is a smoothing coefficient, a very small positive number, and its function is to prevent the denominator from being equal to zero in extreme cases (such as and ) and causing calculation errors, .

[0057] The parameter The acquisition steps are as follows: This parameter represents the total number of non-road feature interference areas identified in the system. After generating the non-road feature interference areas, perform connected component analysis on them, and the total number of independent connected regions counted is , for example, a total of independent interference area patches are identified.

[0058] Calculation process: Take the calculation of the morphological difference interference score of the th candidate unpaved road segment as an example. For example, its area , and through spatial topological indexing or overlap checking, it is found that it only has significant interaction (overlap area greater than 0) with and two interference areas, and the total number of interference areas . The parameters related to the interaction are known as follows: For : , , , . For : For example , it is calculated that , and it is calculated that .

[0059] Calculate the interference term for : ; Calculate the interference term for : ; Calculate the sum (for example, only these two interference areas overlap with ): ; Calculate : ; This result indicates that the morphological difference interference score of the th candidate unpaved road segment is . This score reflects the comprehensive degree of interference of this road segment by the identified non-road features (such as vegetation, water bodies, etc.). The higher the value, the relatively larger the overlap area between this road segment and the interference area, or although the overlap is not large, the shape difference is significant and the distance is close, which increases the possibility that this candidate segment is mis-extracted (for example, the edge of the vegetation adjacent to the road is extracted).

[0060] Based on the morphological difference interference scores of each candidate unpaved road segment calculated in the previous step, screen all candidate road segments, and set an interference score threshold , this threshold needs to be determined through experiments and is aimed at distinguishing real road segments from areas that are mis-extracted due to severe confusion with non-road features. For example, the distribution of known true and false road segment samples can be analyzed, and a boundary that can effectively eliminate most false positives (high values) while retaining most true positives (low values) can be selected. For example, set . Traverse all candidate road segments. If their interference score is greater than (for example ), then remove this road segment from the candidate set and retain all road segments. Next, connect and repair the remaining road segments after screening to handle possible breaks (such as caused by bridges, occlusions, or algorithm limitations). Using the buffer analysis method, first extract the endpoint coordinates of each remaining road segment, and then create a circular buffer with a radius of for each endpoint. This radius is set empirically and is used to define the maximum tolerance gap for endpoint connection. For example, set pixels. Check if there is an overlap between the 5-pixel buffer of one endpoint and the 5-pixel buffer of the endpoint of another different road segment. If there is an overlap, it is considered that these two endpoints represent a break point. A straight connection segment can be added between these two endpoints, or more precisely, the polygons of the two road segments can be merged or geometrically connected. In this way, all break endpoints that meet the conditions are connected, and finally a road network with better connectivity and integrity is formed.

Claims

1. An intelligent road area extraction system for drone aerial images, characterized in that: The system comprises: The image noise filtering module analyzes the frequency domain information of the input drone aerial image, identifies the target noise pattern, filters the drone aerial image, and obtains the filtered image data; The spectral discreteness calculation module divides the drone aerial image into a plurality of local area blocks based on the filtered image data, and calculates the local pixel spectral discreteness value of the pixel color value in each area block; compares and judges the local pixel spectral discreteness value of each area block with a preset discreteness threshold, marks the area below the preset discreteness threshold, and establishes a low spectral discreteness area map; A road connectivity analysis module, based on the low-spectral discrete region map, identifies connected components, obtains connected region morphological parameters, merges low-spectral discrete components according to the connected region morphological parameters, and generates candidate unpaved road segments; The road area constraint refining module identifies vegetation coverage areas and water areas based on the filtered image data, establishes non-road feature interference areas, spatially compares the candidate unpaved road segments with the non-road feature interference areas, removes road segments overlapping with the non-road feature interference areas, and connects the remaining road segments to obtain a road network.

2. The intelligent road area extraction system for drone aerial images according to claim 1 is characterized in that: The steps for obtaining the filtered image data are as follows: The pixel matrix of the input drone aerial image is converted into a frequency domain complex matrix through a two-dimensional fast Fourier transform, and the real component, imaginary component and amplitude spectrum component of the frequency domain complex matrix are extracted to generate frequency domain information including amplitude spectrum kurtosis, phase spectrum variance and frequency energy distribution density; Based on the frequency domain information, traverse the abnormal peak area in the amplitude spectrum component that exceeds 3 times the standard deviation of the average amplitude value, analyze the coordinate position, half-maximum full width and amplitude attenuation gradient of the abnormal peak area, filter the spectrum feature cluster corresponding to the noise according to the coordinate position, and identify the target noise pattern; According to the target noise pattern, a notch filter group is set at the coordinates of the abnormal peak of the amplitude spectrum component, and the center frequency of each notch filter is configured to be the horizontal and vertical coordinate values ​​of the abnormal peak, the bandwidth parameter is 1.5 times the half-maximum full width value, and the attenuation intensity is the inverse of the amplitude attenuation gradient. After filtering the frequency domain complex matrix, an inverse Fourier transform is performed to obtain filtered image data.

3. The intelligent road area extraction system for drone aerial images according to claim 1 is characterized in that: The step of obtaining the local pixel spectrum discreteness value is: The filtered image data is divided into grids with 16×16 pixels as the basic unit, and the sliding step of adjacent grids is set to 8 pixels to generate a set of local area blocks; Based on the local area block set, extract the values ​​of the three channels R, G and B of all pixels in each area block in the RGB color space, calculate the mean, standard deviation, maximum and minimum values ​​of each channel, and generate a color mean vector, a color standard deviation vector, a color maximum vector and a color minimum vector of each area block; The local pixel spectrum dispersion value is calculated based on the color standard deviation vector, the color maximum value vector and the color minimum value vector.

4. The intelligent road area extraction system for drone aerial images according to claim 1 is characterized in that: The steps of obtaining the low-spectrum discrete area map are: Traversing the local pixel spectrum discreteness value sets of all local area blocks, comparing the local pixel spectrum discreteness value of each area block with the preset discreteness threshold one by one, determining the area blocks whose local pixel spectrum discreteness value is less than the preset discreteness threshold, and generating a preliminary marked area index set; Based on the preliminary marked area index set, spatial connectivity analysis is performed on adjacent marked areas, marked areas with a boundary distance less than 2 pixels are merged, and isolated marked areas with an area less than 10 pixels are eliminated to generate an optimized marked area index set; Based on the optimized marked area index set, the coordinates of the marked area are mapped to the original image pixel matrix, the marked area and the unmarked area are filled, and a low-spectrum discrete area map is generated.

5. The intelligent road area extraction system for drone aerial images according to claim 1 is characterized in that: The steps for obtaining the morphological parameters of the connected region are: Based on the low-spectral discrete region map, an eight-neighborhood traversal algorithm is used to traverse the pixels of the entire map, spatially continuous marked regions are aggregated into independent connected components, and a pixel coordinate set and centroid coordinates of each connected component are extracted to generate a connected component set; Based on the connected component set, the Euclidean distances from all pixel coordinates in each connected component to the centroid coordinates are calculated, pixels whose distances are less than 2 times the average Euclidean distance of the components are selected to form a valid direction analysis area, and a set of valid direction analysis areas is generated; Based on the set of valid direction analysis regions, the direction consistency of each connected component is calculated.

6. The intelligent road area extraction system for drone aerial images according to claim 1 is characterized in that: The steps for obtaining the candidate unpaved road segment are as follows: Traversing the morphological parameters of the connected region, screening low-spectral discrete components according to directional consistency, and screening low-spectral discrete components with an aspect ratio greater than 3:1, determining connected components that meet the linear distribution trend, and generating a set of candidate linear components; Based on the candidate linear component set, the distance between the endpoint coordinates of the spatially adjacent components is calculated. If the distance between the adjacent endpoints is less than 2 pixels and the direction angle is less than 15 degrees, the merged component is a continuous area, and an optimized merged component set is generated; Based on the optimized merging component set, the geometric center line of the merging area is extracted, and the boundary of the coverage area is extended along the line direction to generate a convex polygon with uniform width, which is marked as a candidate unpaved road segment.

7. The intelligent road area extraction system for drone aerial images according to claim 1 is characterized in that: The steps for obtaining the non-road feature interference area are as follows: Based on the filtered image data, the reflectance values ​​of each pixel in the near-infrared band and the red band are extracted, the ratio of the near-infrared band reflectance to the red band reflectance is calculated as a preliminary vegetation index, and a preliminary vegetation index grid map is generated; Based on the preliminary vegetation index grid map, the distribution characteristics of the vegetation index of the whole map are counted, and pixels with a vegetation index greater than 0.6 and a near-infrared reflectivity less than 0.3 are set as candidate points of the vegetation coverage area, and pixels with a reflectivity lower than 0.1 and a standard deviation of the blue light band reflectivity less than 0.05 are set as candidate points of the water body area, and a set of vegetation coverage candidate points and a set of water body candidate points are generated; Based on the vegetation coverage candidate point set and the water body candidate point set, the morphological closing operation is used to merge adjacent candidate point areas, isolated areas with an area less than 50 pixels are eliminated, the boundaries of vegetation and water body coverage areas are fused to form continuous patches, and non-road feature interference areas are generated.

8. The intelligent road area extraction system for unmanned aerial vehicle aerial images according to claim 1 is characterized in that: The steps of obtaining the road network are: Based on the candidate unpaved road segment and the non-road feature interference area, the coordinates of the polygonal boundaries of the two are extracted, the spatial topological relationship between the road segment and the interference area is generated by using Delaunay triangulation, and a spatial topological index set is generated; Based on the spatial topological index set, calculating the morphological difference interference score between the candidate road segment and the non-road interference area; Road segments were screened and removed based on the morphological difference interference score, and a 5-pixel buffer was used to analyze the broken endpoints connecting the remaining road segments to generate a road network.

Images