An intelligent extraction system for road areas in UAV aerial images

Abstract

The present invention relates to the field of image processing technology, and specifically to an intelligent road area extraction system for UAV aerial images. The system includes an image noise filtering module that analyzes the frequency domain information of the input UAV aerial image, identifies the target noise pattern, and filters the UAV aerial image. In the present invention, by combining frequency domain noise analysis with 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 the fracture misjudgment under terrains with variable morphologies; 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.

Description

Technical Field

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

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

[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 occlusion or temporary dirt roads, it is difficult to achieve cross-occlusion connection of broken roads based on local features, and a large amount of post-processing interpolation is required, with limited efficiency and accuracy. Therefore, improvements are needed. Summary of the Invention

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

[0005] To achieve the above purpose, the present invention adopts the following technical solutions: An intelligent extraction system for road regions in drone aerial images includes:

[0006] An image noise filtering module, which 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;

[0007] A spectral dispersion calculation module, based on the filtered image data, divides the drone aerial image into multiple local region blocks, calculates the local pixel spectral dispersion values 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;

[0008] A road connectivity analysis module, based on the low spectral dispersion region map, identifies the 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;

[0009] The road area constraint refinement module, based on the filtered image data, identifies the vegetation-covered area and the water area, establishes a non-road feature interference area, compares the spatial positions of the candidate unpaved road segments with the non-road feature interference area, removes the road segments that overlap with the non-road feature interference area, and connects the remaining road segments to obtain a road network.

[0010] Preferably, the steps for obtaining the filtered image data are as follows:

[0011] Convert the pixel matrix of the input UAV aerial image into a frequency-domain complex matrix through two-dimensional fast Fourier transform, extract the real component, imaginary 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;

[0012] 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;

[0013] 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 to be the horizontal and vertical coordinate values of the abnormal peak, the bandwidth parameter to be 1.5 times the full width at half maximum value, and the attenuation intensity to be the reciprocal of the amplitude attenuation gradient. After filtering the frequency-domain complex matrix, perform inverse Fourier transform to obtain the filtered image data.

[0014] Preferably, the steps for obtaining the local pixel spectral dispersion value are as follows:

[0015] Divide the filtered image data into grids with 16×16 pixels as the basic unit, set the sliding step of adjacent grids to be 8 pixels, and generate a set of local region blocks;

[0016] Based on the set of local region blocks, extract the numerical values of the R, G, and B channels of all pixels in each region block in the RGB color space, calculate the mean value, standard deviation, maximum value, and minimum value of each channel, and generate the color mean vector, color standard deviation vector, color maximum vector, and color minimum vector of each region block;

[0017] Based on the color standard deviation vector, color maximum vector, and color minimum vector, calculate the local pixel spectral dispersion value.

[0018] Preferably, the steps for obtaining the low spectral dispersion area map are as follows:

[0019] Traverse the set of local pixel spectral dispersion values of all local region blocks, compare each local pixel spectral dispersion value of the region block item by item with a preset dispersion threshold, determine the region blocks where the local pixel spectral dispersion value is less than the preset dispersion threshold, and generate a preliminary marked region index set;

[0020] Based on the preliminary marked region index set, perform spatial connectivity analysis on adjacent marked regions, merge the marked regions with a boundary distance less than 2 pixels, and eliminate the isolated marked regions with an area less than 10 pixels to generate an optimized marked region index set;

[0021] Based on the optimized marked region index set, map the coordinates of the marked regions to the original image pixel matrix, fill the marked regions and non-marked regions, and generate a low spectral dispersion region map.

[0022] Preferably, the steps for obtaining the morphological parameters of the connected regions are as follows:

[0023] Based on the low spectral dispersion region map, traverse all the pixels of the image using the eight-neighborhood traversal algorithm, aggregate the 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;

[0024] Based on the connected component set, calculate the Euclidean distance from all pixel coordinates in each connected component to the centroid coordinate, screen the pixels with a distance less than 2 times the average Euclidean distance of the component to form an effective direction analysis region, and generate an effective direction analysis region set;

[0025] Based on the effective direction analysis region set, calculate the direction consistency of each connected component.

[0026] Preferably, the steps for obtaining the candidate non-paved road segments are as follows:

[0027] Traverse the morphological parameters of the connected regions, screen the low spectral dispersion components according to the direction consistency, and at the same time screen the low spectral dispersion components with an aspect ratio greater than 3:1, determine the connected components that conform to the linear distribution trend, and generate a candidate linear component set;

[0028] Based on the candidate linear component set, calculate the distance 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;

[0029] Based on the optimized merged component set, extract the geometric center connection line of the merged region, extend 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.

[0030] Preferably, the steps for obtaining the non-road ground object interference area are as follows:

[0031] 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;

[0032] Based on the preliminary vegetation index raster map, statistically analyze the distribution characteristics of the vegetation index across the entire 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;

[0033] 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 regions, remove isolated regions with an area less than 50 pixels, and fuse the boundaries of the vegetation and water-covered areas to form continuous patches, generating a non-road ground object interference area.

[0034] Preferably, the steps for obtaining the road network are as follows:

[0035] Based on the candidate unpaved road segments and the non-road ground object interference area, 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 area, generating a set of spatial topological indexes;

[0036] Based on the set of spatial topological indexes, calculate the morphological difference interference score between the candidate road segments and the non-road interference area;

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

[0038] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0039] In the present invention, by combining frequency domain noise analysis with notch filtering, high-frequency noise and periodic interference are suppressed while retaining the sharpness of the 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, combined with direction consistency and aspect ratio constraints, enhancing the connectivity perception ability of fragmented unpaved roads and reducing the misjudgment of fractures under terrains with diverse morphologies; using the normalized vegetation index and multi-band reflectance for joint modeling, excluding water body reflection interference through the standard deviation of the blue band, and constructing a spatial mask for the non-road ground object interference area; based on the morphological difference scoring mechanism of the Hausdorff distance and shape similarity index, quantifying the multi-dimensional spatial relationship between the road segments and the interference area, and retaining candidate road segments with independent morphology and slight contact when removing overlapping areas, reducing the over-trimming problem caused by the single index of the traditional area overlap rate. Brief Description of the Drawings

[0040] Figure 1 This is the system flowchart of the present invention. Detailed Description of the Invention

[0041] In order to make the objectives, technical solutions and advantages of the present invention clearer and more understandable, the present invention will be further described in detail below with reference to the 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.

[0042] Please refer to Figure 1 , the present invention provides a technical solution: An intelligent extraction system for road areas in drone aerial images includes:

[0043] An 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;

[0044] A spectral dispersion calculation module divides the drone aerial image into multiple local region blocks based on the filtered image data, calculates the local pixel spectral dispersion values of the pixel color values within each region block; compares the local pixel spectral dispersion values of each region block with a preset dispersion threshold for judgment, marks the regions below the preset dispersion threshold, and establishes a low spectral dispersion region map;

[0045] A road connectivity analysis module identifies connected components based on the low spectral dispersion region map, 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;

[0046] A 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.

[0047] The steps for obtaining the filtered image data are as follows:

[0048] Convert the pixel matrix of the input drone 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;

[0049] Based on the frequency-domain information, traverse the abnormal peak regions in the amplitude spectrum components that exceed 3 times the standard deviation of the average amplitude value, analyze the coordinate positions, full width at half maximum (FWHM), and amplitude attenuation gradient of the abnormal peak regions, screen the spectral feature clusters corresponding to noise according to the coordinate positions, and identify the target noise patterns;

[0050] According to the target noise patterns, set a notch filter bank at the abnormal peak coordinates of the amplitude spectrum components. 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 FWHM value, and the attenuation intensity as the reciprocal of the amplitude attenuation gradient. After filtering the frequency-domain complex matrix, perform the inverse Fourier transform to obtain the filtered image data.

[0051] 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, generating 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 respectively to form two independent real matrices, namely the real component matrix and the imaginary 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 , thereby 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 , and its variance reflects the consistency or dispersion of the phase information. Finally, according to 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 sum or average of the energy in the low-frequency, medium-frequency, and high-frequency regions (which can be divided according to the frequency radius. For example, divide the frequency-domain origin radius 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 amplitude spectrum kurtosis, phase spectrum variance, and frequency energy distribution density to obtain the frequency-domain information characterizing the image noise characteristics.

[0052] 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 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 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 within 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 i.e., the frequency width corresponding to when the amplitude drops to half of the peak, and estimate the attenuation gradient of the amplitude value near the peak center with respect to frequency, for example, obtained by calculating the average gradient of the amplitude values within the neighborhood of the peak center. Finally, based on the coordinate positions of these extracted abnormal peak regions, 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.

[0053] According to the specific information (coordinates, full width at half maximum, amplitude attenuation gradient) of each abnormal peak included 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 notch filter bank. When configuring specifically, set the center frequency of each notch filter to the center coordinate value of the corresponding abnormal peak. The bandwidth parameter of the filter is set according to the full width at half maximum (FWHM) of this peak, specifically set to 1.5 times the measured full width at half maximum. For example, if the full width at half maximum of a certain 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 the direction by units, and the attenuation intensity of the filter is determined according to the reciprocal of the amplitude attenuation gradient of this 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, and all configured individual notch filters will be combined into a total filter transfer function , the value of this function is less than 1 in the notch region and 1 in other regions. Subsequently, this total filter transfer function is multiplied element by element with the original frequency-domain complex matrix to obtain the filtered frequency-domain complex matrix. Finally, a two-dimensional inverse fast Fourier transform is performed on this filtered frequency-domain complex matrix to convert it from the frequency domain back to the spatial domain, obtaining the filtered image data.

[0054] The steps to obtain the local pixel spectral dispersion value are as follows:

[0055] Divide the filtered image data into grids with a basic unit of 16×16 pixels, set the sliding step of adjacent grids to 8 pixels, and generate a set of local region blocks;

[0056] 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;

[0057] Based on the color standard deviation vector, color maximum value vector, and color minimum value vector, calculate the local pixel spectral dispersion value. The calculation formula is: ;

[0058] where is the local pixel spectral dispersion value of the s-th region 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 means of the other two channels.

[0059] 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 multiplied by 16 pixels as the size of each local region block, traverse the entire filtered image data, start from the top-left coordinate (0, 0) of the image, extract the first 16×16 pixel block, and then set the sliding steps in both the horizontal and vertical directions to 8 pixels. This means that the next horizontally adjacent block will start extracting 16×16 pixels from the coordinate (8, 0), 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 Less than the image width), similarly, after completing the extraction of the first row of blocks, move to the next row. The starting block coordinates are (0, 8), and the subsequent blocks are (8, 8), (16, 8), etc., until the entire image area is covered (the starting row coordinate of the last block Meet Less than the image height). This sliding window method with overlap ensures a smooth transition of pixel neighborhood information between blocks. Organize the coordinates of all the extracted 16×16 pixel blocks and their corresponding pixel data to generate a set of local area blocks.

[0060] Based on the set of local area blocks generated in the previous step, perform color feature statistics on each local area block in the set (i.e., a 16×16 pixel block). First, access all 256 pixel points within the area block and read the three-channel values of each pixel point in the RGB color space, namely the red (R) channel value, the green (G) channel value, and the blue (B) channel value. These values are usually in the integer range of 0 to 255. Then, for the 256 pixels within the area 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 value to obtain the average color of the block. Calculate the standard deviation of each channel's numerical value 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 value. Subsequently, combine the calculated average values of the three channels 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 (for example ), combine the maximum values of the three channels into a color maximum vector (for example ), and combine the minimum values of the three channels into a color minimum vector (for example ). Repeat this process for all local area blocks, and finally generate the color mean vector, color standard deviation vector, color maximum vector, and color minimum vector corresponding to each area block.

[0061] Formula: , The benefit of the formula is that it comprehensively considers the changes within each color channel (standard deviation and range ) and the mutual relationship between channels (covariance ), so as to more accurately quantify the spectral uniformity or complexity of the local area. 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 mean values of the other two channels. When the color channels in a region change in concert (for example, when the overall brightness rises or falls, causing the R, G, and B values to increase or decrease synchronously, and the absolute value of the covariance is large), it indicates that the 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. On the contrary, if the change relationship between channels is chaotic (the absolute value of the covariance is small), the greater influence of the range is retained, which can better reflect the characteristics of mixed pixels or complex texture regions. This design enables to give lower values to areas such as roads with relatively single spectra and relatively stable internal color relationships, and higher values to complex areas such as vegetation and building edges, improving the discrimination;

[0062] Parameter is obtained as follows: This parameter represents the standard deviation of the pixel values of the color channel within the ( is one of R, G, B), which 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, it is directly extracted from its color standard deviation vector . If the vector value is , then , , .

[0063] Parameter is obtained as follows: This parameter represents the range of the pixel values of the color channel within the th region block, that is, the difference between the maximum value and the minimum value of this channel, which is calculated from the color maximum value vector and the color minimum value vector . 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 .

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

[0065] Calculation process: Taking the th local region block as an example, substitute the parameter values obtained previously for calculation: Known parameters: , , , , , , ;

[0066] Calculate the contribution of the R channel: ;

[0067] Calculate the contribution of the G channel: ;

[0068] Calculate the contribution of the B channel: ;

[0069] Calculate the total local pixel spectral dispersion value : ;

[0070] This result indicates that the calculated result of the local pixel spectral dispersion value of the th local region block is . 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, which needs to be determined based on the distribution statistics of road and non-road regions in a large number of sample images, and select a boundary value that can better distinguish the two, such as setting it to 30) tends to indicate that the spectral characteristics of this region are relatively single and stable, and may be part of homogeneous regions such as roads and bare lands. A higher value (for example, greater than 30) indicates that the spectral complexity of this region is high, and may be a reflection of vegetation, buildings, mixed landforms, etc.

[0071] The steps to obtain the low spectral dispersion region map are as follows:

[0072] Traverse the set of local pixel spectral dispersion values of all local region blocks, compare each local pixel spectral dispersion value of each region block with a preset dispersion threshold item by item, determine the region blocks with local pixel spectral dispersion values less than the preset dispersion threshold, and generate a preliminary marked region index set;

[0073] Based on the preliminary marked region index set, perform spatial connectivity analysis on adjacent marked regions, merge the marked regions with a boundary distance less than 2 pixels, and remove the isolated marked regions with an area less than 10 pixels to generate an optimized marked region index set;

[0074] Based on the optimized marked region index set, map the coordinates of the marked regions to the original image pixel matrix, fill the marked regions and unmarked regions to generate a low spectral dispersion region map.

[0075] Specifically, traverse the set of local pixel spectral dispersion values of all local region blocks obtained in the previous step The set contains the dispersion values calculated for each 16×16 pixel region block. For each region block value is compared with a preset dispersion threshold The determination process of this preset dispersion threshold is as follows: Select a representative set of UAV aerial image samples covering various scenarios (different lighting, weather, surface cover). Manually and precisely outline typical road regions (especially unpaved roads) and non-road regions (such as vegetation, buildings, water bodies, bare land, etc.) in the sample images. Apply the same 16×16 overlapping block division and calculation process to these manually labeled regions to obtain a large number of values of known categories values. Statistically analyze the distribution histograms of the values in the road regions values and non-road regions values, analyze the overlap situation of the two distributions, and select a value that can better distinguish the two types of regions value as the threshold. For example, the Otsu method can be used to automatically calculate a threshold, or based on empirical observation, select a value located between the peaks of the two distributions. For example, if statistical analysis shows that most of the values of the road region blocks are concentrated between 15 and 35, while most of the values of the non-road region blocks are distributed above 40, then the preset dispersion threshold can be set to 30. During the traversal process, for the th region block, if its value is less than value less than For example ), it is determined that the regional block belongs to the low-spectral discrete region, and the index of the regional block (such as its row, column number or starting pixel coordinates in the grid) is recorded. The indices of all regional blocks that meet the conditions are gathered to generate a preliminary marked region index set.

[0076] Based on the preliminary marked region index set generated in the previous step, optimize the spatial relationship of these regional blocks preliminarily marked as low-spectral discrete. First, map the position information of these marked regional blocks to a two-dimensional space (corresponding to the original image space), regard each marked block as a spatial entity, conduct spatial connectivity analysis, and find spatially adjacent or close marked regional blocks. 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 empirically 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 merge these two regions 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 (by de-duplicating and counting the pixel points contained in all 16x16 blocks that make up the region), and set a minimum area threshold , which is used to remove low-discrete regions that are too small and are likely to be noise or accidentally formed by non-road features. This threshold is determined empirically based on the minimum effective width and length of the road and the image resolution in practical applications. For example, for an application aimed at extracting the carriageway, it 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 it is removed from the marked regions. Retain all connected regions with an area greater than or equal to . The corresponding regional block indices form an optimized marked region index set.

[0077] Based on the optimized marked region index set generated in the previous step, which contains the region block indexes that are considered 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 the same size as the original filtered image data (e.g., 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 the 16×16 region block in the original image coordinate system according to the index (e.g., the block index corresponds to the upper-left pixel coordinate , and the coverage range is ). Next, modify the pixel values at the corresponding positions of all the pixels covered by this region block (a total of 256) 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 pixels 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 the 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, obtaining the low-spectral discrete region map.

[0078] The steps for obtaining the morphological parameters of connected regions are as follows:

[0079] Based on the low-spectral discrete region map, traverse all the pixels of the whole image using the eight-neighborhood traversal algorithm, aggregate the 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;

[0080] 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, generating an effective direction analysis region set;

[0081] Based on the effective direction analysis region set, calculate the direction consistency of each connected component. The calculation formula is: ;

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

[0083] 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), identify and separate each independent marked region in the map. Specifically, use the eight-neighborhood traversal algorithm. This algorithm scans row by row starting from the first pixel of the image (e.g., the upper left corner). When encountering an unvisited marked pixel (with a value of 1), it is used as the starting point (seed point) of a new connected component and is assigned a unique identifier for the component. At the same time, initialize the pixel coordinate list, coordinate sum (for calculating the centroid), and pixel counter of the component. Then, use a queue or stack structure to perform breadth-first search or depth-first search. Add the seed point to the queue / stack and mark it as visited. Loop through the pixels in the queue / stack: Take out a pixel and 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, mark it as visited, assign it the same component identifier, add its coordinates to the pixel coordinate list of the current component, accumulate its coordinate values to the sum, increase the pixel counter, and add it 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 coordinate sum 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 have been visited and belong to a connected component. Finally, generate a set of connected components containing information about all independent connected components.

[0084] 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 better represents 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 , then, 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 eliminating 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 as 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 as a valid pixel. All the retained pixels together constitute the effective direction analysis area of the th connected component. Perform the same operation on all connected components to generate a set of effective direction analysis areas.

[0085] Formula: , the benefit of the formula is that it provides a method to measure the consistency of the distribution direction of the 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 extending bidirectionally along a single axis (such as a straight line segment);

[0086] The steps to obtain the parameter are as follows: This parameter represents the total number of pixels contained 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 .

[0087] 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 The centroid of a component is , and the th pixel coordinate in its effective direction analysis region is , then the azimuth angle of this pixel is radians.

[0088] 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 th connected component. These coordinates are directly derived 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 .

[0089] 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 , .

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

[0091] Calculation : ;

[0092] Calculate : ;

[0093] Calculate : ;

[0094] The result shows that for this simplified example connected component with 3 valid pixels , its direction consistency The calculated value is , and this value is between 0 and 1. The closer the value is to 1, it indicates that the pixel points forming the effective area tend to be distributed on a straight line passing through the centroid, showing a strong linear feature, which is usually the morphological manifestation of slender ground objects such as roads. The closer the value is to 0, it indicates that the pixel points are more diffusely distributed, lacking obvious directionality, and are more like blocky areas.

[0095] The steps to obtain candidate unpaved road segments are as follows:

[0096] 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, determine the connected components that conform to the linear distribution trend, and generate a set of candidate linear components;

[0097] Based on the set of candidate linear components, calculate the distance 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;

[0098] Based on the optimized merged component set, extract the geometric center connection line of the merged region, extend the boundary of the coverage area along the connection line direction to generate a convex polygon with a uniform width, and mark it as a candidate unpaved road segment.

[0099] Specifically, traverse each connected component obtained in the previous steps and its associated morphological parameters, especially the direction consistency value, filter these low-spectral discrete components, and set a direction consistency threshold , which is determined according to prior knowledge or statistical analysis of sample data, aiming to retain the components with obvious linear features. For example, by analyzing the value distribution of known road samples and non-road sample components, it is found that the value of road components is usually higher than 0.7, while the value of non-road components (such as blocky farmland, parking lots) is lower. Therefore, it can be set , if the The value is greater than (for example ), it is initially considered to meet the requirements of linear characteristics. At the same time, calculate the aspect ratio of the component. First, use the principal component analysis method to calculate the major axis direction and minor axis direction of the set of pixel coordinates of the component, and then calculate the maximum projection length of the component in the major axis direction 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, this threshold is set to 3:1, that is . If the aspect ratio of the component is greater than (for example, calculated as ), it is considered to conform to the slender characteristic in shape. Only the connected components that simultaneously meet the two conditions that the direction consistency is greater than and the aspect ratio is greater than are determined to conform to the linear distribution trend and are selected and added to the candidate linear component set

[0100] 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 major 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 (for example, 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 major 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), and calculate the Euclidean distance between the closest pair of endpoints between them , and at the same time calculate the included angle (take the acute or obtuse angle between 0 and 180 degrees) between the main direction vectors of these two components. Set two merge condition thresholds: the endpoint distance threshold and the direction included angle threshold . These two thresholds are set according to experience to allow connecting broken road segments with a small gap 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 included angle is less than (for example ), it is determined that these two components should be merged. The pixel sets of them are merged, and the relevant attributes of the merged region are updated (such as recalculating the end points and the main direction). Repeat this process, using, for example, graph-based connection methods, until all adjacent components that meet the conditions are merged, generating an optimized merged component set.

[0101] Based on the optimized merged component set generated in the previous step, where each element represents a potentially more continuous road area after merging, perform geometric normalization on each merged region. First, extract the connecting line of the geometric centers of the merged region, 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 region, or skeletonizing the region to extract the longest skeleton path, or simply connecting the centroids of the original (before merging) candidate linear components included in the merged region and smoothing to obtain a polyline. 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 perpendicular directions on its two sides by distance (that is, extend 5 pixels on each side), forming a series of equal-width line segments perpendicular to the center line. All the end points of these line segments outline the set of boundary contour points of the road region. Finally, apply the convex hull algorithm to this set of boundary contour points, 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.

[0102] The steps to obtain the non-road feature interference area are as follows:

[0103] 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;

[0104] 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 in 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 in the water area, generating a set of candidate points in the vegetation-covered area and a set of candidate points in the water area;

[0105] Based on the vegetation coverage candidate point set and the water body candidate point set, morphological closing operation is used to merge adjacent candidate point regions, isolated regions with an area less than 50 pixels are removed, the boundaries of the vegetation and water body coverage areas are fused to form continuous patches, and a non-road ground object interference area is generated.

[0106] Specifically, based on, for example, available filtered image data containing near-infrared and red bands (the filtered image data here refers to multi-spectral or hyperspectral image data that has undergone preprocessing such as radiometric calibration and atmospheric correction to obtain surface reflectance values), an operation is performed on each pixel in the image. First, for pixel position , the reflectance value in the near-infrared band is extracted respectively and the reflectance value in the red light band . These reflectance values are dimensionless and usually range from 0 to 1. Subsequently, the ratio of these two reflectance values is calculated, 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 maximum value or invalid value mark. This ratio is used as a simplified preliminary vegetation index, and the magnitude of its value is usually positively correlated with the density of surface vegetation coverage. Because healthy vegetation strongly reflects near-infrared light and strongly absorbs red light, the values calculated for all pixels in the image are stored in a new two-dimensional raster data structure, which has the same spatial resolution and range as the original image, generating a preliminary vegetation index raster map.

[0107] Based on the preliminary vegetation index raster map generated in the previous step and combined with the original filtered image data containing reflectance in bands such as near-infrared, red, and blue, pixels are classified to identify vegetation and water body candidate regions. First, the distribution characteristics of all pixel values in the preliminary vegetation index raster map are statistically analyzed, such as calculating its histogram, mean, and standard deviation, to understand the approximate range of the vegetation index in the study area. The determination criteria for candidate points in the vegetation coverage area are set: 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 some building roofs) or specific types of vegetation (such as very dense and dark-colored vegetation). For candidate points in the water body area, the determination criteria are set: its reflectance needs to be lower than the threshold , and the standard deviation of the reflectance in the blue light band within the local neighborhood needs to be less than the threshold , where the reflectance in the blue light band can be selected , because the reflectance of water bodies is usually the lowest in this band, and the threshold is set to 0.1 based on experience, reflecting the low reflectance characteristics of water bodies; the local neighborhood standard deviation is obtained by calculating the standard deviation of the reflectance in the blue light band of all pixels within a window such as a 5×5 window around each pixel. This value is used to measure the texture uniformity of the region. The water surface is usually very uniform, so is set to a small value, such as 0.05, to distinguish it from shadows or asphalt roads with the same low reflectance but possible textures. Traverse each pixel in the image and make judgments according to the above criteria: If the vegetation conditions are met (such as and ), then record its coordinates as vegetation coverage candidate points; if the water body conditions are met (such as and ), then record its coordinates as water body area candidate points, and finally generate a set of vegetation coverage candidate points and a set of water body candidate points.

[0108] 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, map the set of vegetation candidate points and the set of water body candidate points to two independent binary images respectively. Set the candidate point positions to the foreground value 1 and the rest to the background value 0. Perform morphological closing operations on these two binary images independently. The closing operation consists of dilation followed by erosion. Select a small structuring element, such as a 3×3 pixel-sized all-1 square kernel. The size of this structuring element is selected according to experience, aiming to fill small gaps between candidate points and small holes inside the region, so that originally close but disconnected candidate point regions can be connected to form more complete patches. After the closing operation, perform connected component analysis (using eight-neighborhood connection) on the two binary images again to identify all independent connected regions composed of foreground pixels, calculate the pixel area of each connected region, and set a minimum area threshold , which is set based on experience and is used to filter out those small-area, likely noisy or insignificant sporadic patches. For example, set For a pixel, if the area of a connected region is less than 50 pixels, then set all pixel values within that region from 1 to 0, and eliminate it. Retain all regions with an area greater than or equal to 50 pixels. Finally, perform a logical OR operation on the binary vegetation region map and the binary water body region map after area filtering, that is, create a new blank binary map. If a pixel is foreground value 1 in at least one of the processed vegetation map or the processed water body map, then set that pixel to 1 in the new map, otherwise set it to 0. In this way, the verified vegetation areas and water body areas are merged into continuous patches to generate the final non-road ground object interference area.

[0109] The steps to obtain the road network are as follows:

[0110] Based on the candidate unpaved road segments and the non-road ground object interference area, 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 area, generating a spatial topological index set;

[0111] Based on the spatial topological index set, calculate the morphological difference interference score between the candidate road segments and the non-road interference area. The calculation formula is: ;

[0112] Where, 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 of the Hausdorff distance; , is the shape similarity index, and are the perimeters of the road segment and the interference area respectively, is the pixel area of the k-th candidate unpaved road segment, is the pixel area of the m-th non-road ground 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 ground object interference areas;

[0113] Filter and remove the road segments according to the morphological difference interference score, and use a 5-pixel buffer analysis to connect the broken endpoints of the remaining road segments to generate the road network.

[0114] Specifically, based on the candidate unpaved road segments and non-road feature 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 feature interference area polygon, merge all these vertex coordinates from the road segments and interference areas into a unified vertex set, and 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 examining the edges of this 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) that have this topological proximity relationship to generate a spatial topological index set.

[0115] Formula: , The benefit of the formula is that it provides a quantitative indicator for evaluating the degree to which each candidate unpaved road segment is affected by the non-road feature interference area , This indicator 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 their boundaries. The farther the distance, even if there is overlap, its 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 (such as a slender road and a circular interference area), even if the distance is close and there is overlap, it may indicate that the interference area has little association with the road structure, and by increasing to reduce its influence weight;

[0116] The steps to obtain the 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 this polygon. For example, after calculation, the th candidate road segment polygon covers 1500 pixels, then .

[0117] Parameter is obtained 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, if the area of the th interference area (such as a piece of vegetation) is calculated to be 6200 pixels, then .

[0118] Parameter is obtained as follows: This parameter represents the overlapping pixel area between the th candidate unpaved road segment and the th non-road feature interference area. The spatial intersection of the two areas (polygon or raster) needs to be calculated. If both are polygons, the intersection polygon can be obtained using polygon clipping or Boolean operations in computational geometry, and then its area is calculated; for example, if it is calculated that the th road segment overlaps with the th interference area by 450 pixels, then .

[0119] Parameter is obtained as follows: This parameter represents the modified Hausdorff distance between the boundary point set of the th road segment and the boundary point set of the th 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 .

[0120] Parameter is obtained as follows: This parameter represents the shape (compactness) similarity index 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 areas and as before, then calculate their perimeter-area ratios and , finally take the absolute value of the difference between these two ratios, i.e., . The smaller this value is, the more similar the compactness of the two shapes. For example, for the th road segment with a perimeter of pixels and an area of pixels; for the th interference area with a perimeter of pixels and an area of pixels, then .

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

[0122] Parameter is obtained 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 is . For example, a total of independent interference area patches are identified.

[0123] Calculation process: Take the calculation of the morphological difference interference score of the th candidate non-paved road segment as an example. For example, its area is , 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 is . The parameters related to the interaction are known as follows: For : , , , . For : For example, , the calculated , and the calculated .

[0124] Calculate the interference term for : ;

[0125] Calculate the interference term for : ;

[0126] Calculate the sum (for example, only these two interference areas are related to Overlap):

[0127] ;

[0128] Calculate :

[0129] ;

[0130] The result shows that the morphological difference interference score of the th candidate unpaved road segment is , and 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 overlapping area of this road segment with the interference area, or although the overlap is not large, the shape difference is significant and the distance is very 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).

[0131] According to the morphological difference interference scores of each candidate unpaved road segment calculated in the previous step, screen all candidate road segments, set an interference score threshold , and this threshold needs to be determined through experiments, aiming to distinguish real road segments from areas mis-extracted due to serious 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) is selected. For example, set . Traverse all candidate road segments. If its 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). Adopt 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 according to experience and is used to define the maximum tolerance gap for endpoint connection. For example, set For each pixel, check if there is an overlap between the 5-pixel buffer of an endpoint and the 5-pixel buffer of an endpoint of a different road segment. If there is an overlap, it is considered that these two endpoints represent a break point. A straight line 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. By this way, all the break endpoints that meet the conditions are connected, and finally a road network with better connectivity and integrity is formed.

Claims

1. An intelligent extraction system for road areas in UAV aerial images, characterized in that, The system includes: An image noise filtering module that 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 that divides the UAV aerial image into multiple local region blocks based on the filtered image data, calculates the local pixel spectral dispersion values of the pixel color values within each region block; compares and determines 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; A road connectivity analysis module that identifies connected components based on the low spectral dispersion region map, 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 region constraint refinement module that identifies vegetation-covered areas and water areas based on the filtered image data, establishes non-road feature interference areas, compares the spatial positions of the candidate unpaved road segments with the non-road feature interference areas, removes the road segments that overlap with the non-road feature interference areas, and connects the remaining road segments to obtain a road network; 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; The steps for obtaining the local pixel spectral dispersion values are as follows: Divide the filtered image data into grids with a basic unit of 16×16 pixels, set the sliding step of adjacent grids as 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 within 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 vector, and color minimum 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 values; The steps for obtaining the low spectral dispersion region map are as follows: Traverse the set of local pixel spectral dispersion values of all local region blocks, compare each local pixel spectral dispersion value of the region block item by item with the preset dispersion threshold, determine the region blocks where the local pixel spectral dispersion value is less than the preset dispersion threshold, and generate a preliminary marked region index set; Based on the preliminary marked region index set, perform spatial connectivity analysis on adjacent marked regions, merge the marked regions with a boundary distance less than 2 pixels, and remove the isolated marked regions with an area less than 10 pixels to generate an optimized marked region index set; Based on the optimized marked region index set, map the coordinates of the marked regions to the original image pixel matrix, fill the marked regions and non-marked regions, and generate a low spectral dispersion region map.

2. The intelligent extraction system for road regions in UAV aerial images according to claim 1, wherein, The steps for obtaining the morphological parameters of the connected regions are as follows: Based on the low spectral dispersion region map, traverse all the image pixels using the eight-neighborhood traversal algorithm, aggregate the 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, filter out the pixels with a distance less than 2 times 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.

3. The intelligent extraction system for road regions in UAV aerial images according to claim 1, wherein 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 dispersion components according to the direction consistency, and at the same time filter out the low spectral dispersion components with an aspect ratio greater than 3:1, determine the connected components that conform to the linear distribution trend, and generate a candidate linear component set; Based on the candidate linear component set, calculate the distance 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, extend 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.

4. The intelligent extraction system for road areas in UAV aerial images according to claim 1, characterized in that, 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 of the whole image. Set the pixels with a vegetation index greater than 0.6 and a near-infrared reflectance less than 0.3 as candidate points in the vegetation coverage 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 in the water body area, and generate a vegetation coverage candidate point set and a water body candidate point set; Based on the vegetation coverage candidate point set and the water body candidate point set, 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.

5. The intelligent extraction system for road areas in UAV aerial images according to claim 1, characterized in that, 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, 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; Screen and remove the road segments according to the morphological difference interference score, and use 5-pixel buffer analysis to connect the broken endpoints of the remaining road segments to generate a road network.

Images