An automatic geographic information extraction method based on remote sensing images

Abstract

The application relates to the technical field of remote sensing image processing and geographic information extraction, and particularly discloses an automatic geographic information extraction method based on remote sensing images, which comprises the following steps: acquiring remote sensing images, rational polynomial parameters and digital elevation data; generating a binary raster; calculating the slope and slope direction according to the digital elevation data, and non-uniformly stretching the target ground object image element in the binary raster along the slope direction in combination with the observation view angle to obtain a compensation raster and a stretching amplitude; tracking the boundary of the compensation raster, extracting contour nodes and synchronously acquiring the stretching amplitude of each node to form a contour node sequence with the stretching attribute; dynamically determining a simplified threshold according to the stretching attribute of each node, eliminating the nodes with the vertical distance less than the threshold, obtaining a simplified contour and converting the simplified contour into a geographic coordinate vector; and the application eliminates the terrain perspective compression, and makes the contour simplification adapt to the terrain fluctuation.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing and geographic information extraction technology, specifically to an automatic geographic information extraction method based on remote sensing images. Background Technology

[0002] With the rapid development of high-resolution remote sensing satellite technology, automatic extraction of geographic information from remote sensing images has become an important technical means in fields such as land surveys, forestry monitoring, and urban planning. Existing automatic geographic information extraction methods typically include steps such as semantic segmentation, boundary tracking, polyline simplification, and vectorization output.

[0003] In existing methods for automatic geographic information extraction based on remote sensing images, the polyline simplification algorithm with a fixed vertical distance threshold cannot distinguish between pixel-level dense fluctuations caused by topographic perspective contraction and the actual minute transitions of features conforming to the terrain when processing the contours of features in mountainous areas. As a result, contour nodes that reflect the topographic undulation characteristics are discarded as redundant points, and the simplified vector contours deviate from the actual spatial positions of features, making it impossible to perform accurate overlay analysis with digital elevation models. Summary of the Invention

[0004] The purpose of this invention is to provide an automatic geographic information extraction method based on remote sensing imagery to solve the problems mentioned above.

[0005] The objective of this invention can be achieved through the following technical solutions: An automatic geographic information extraction method based on remote sensing imagery includes the following steps: S1, acquire remote sensing images of the target area, rational polynomial parameters corresponding to the remote sensing images, and digital elevation data of the target area; S2 performs pixel-by-pixel feature attribution determination on remote sensing images and generates a binary raster that records the spatial distribution of target features. S3. Based on digital elevation data, the terrain slope and aspect are calculated pixel by pixel. Combined with the observation angle in the rational polynomial parameters, the pixels belonging to the target features in the binary grid are non-uniformly stretched along the aspect direction to restore the width of the compressed features to the horizontal projection width, thus obtaining a compensation grid and recording the stretching amplitude of each pixel to form a compensation data set. S4, perform boundary tracking on the compensation grid in the compensation data group, extract the initial contour nodes of the target features, and simultaneously obtain the stretching amplitude of each node in the compensation data group to form a sequence of contour nodes with stretching attributes. S5 dynamically determines the simplification threshold of each node based on its own stretching properties in the contour node sequence, removes nodes with vertical distances less than the simplification threshold, obtains the simplified contour, and converts the simplified contour into a geographic coordinate vector.

[0006] As a further aspect of the present invention: the generation of a binary raster recording the spatial distribution of target features specifically includes: Select the values ​​of the target ground features in at least two spectral bands in the remote sensing image, calculate the value difference of the same pixel in different spectral bands, and arrange the obtained differences according to the pixel position to form the first intermediate grid indicating the degree of difference between the target ground features and other ground features in spectral response. For each pixel in the first intermediate grid, the values ​​of all pixels in its neighborhood are taken to form a local numerical sequence. The dispersion value of the local numerical sequence is calculated and compared with a preset dispersion threshold. Based on the comparison result, the value of the corresponding pixel in the first intermediate grid is corrected to weaken isolated noise pixels and retain continuous edge regions, thus forming the second intermediate grid. The value of each pixel in the second intermediate grid is compared with a preset binarization threshold. Pixels that are greater than or equal to the preset binarization threshold are marked as target features, and pixels that are less than the preset binarization threshold are marked as background, thus generating a binary grid.

[0007] As a further aspect of the present invention: the composition of the compensation data group specifically includes: The maximum rate of change of the elevation difference between each pixel and its neighboring pixels in the digital elevation data is calculated in both the horizontal and vertical directions to determine the slope and aspect of the pixel. The slope aspect and the observation angle are compared in space. The spatial position relationship between the slope surface where the pixel is located and the observation direction of the sensor is determined based on the angle between the two. Then, the main direction of the pixel affected by perspective shrinkage is determined, and the stretching baseline is set along the main direction. Based on the slope, the number of raster positions that need to be expanded for each cell is determined in the vertical direction of the stretching baseline. The attribute values ​​of the cells are then migrated to the new raster positions along the stretching baseline according to the number of raster positions. At the same time, the actual number of raster positions migrated for each cell is saved as the stretching amplitude of the cell. After traversing all target ground objects, a compensation raster and compensation data set are generated.

[0008] As a further aspect of the present invention: determining the principal direction of the pixel affected by perspective shrinkage and setting a stretching baseline along the principal direction specifically includes: Using the location of the pixel as a reference, the angle values ​​of the slope aspect and the angle values ​​of the observed line of sight are mapped to the same circular measurement system, and the angle between the two in the circular measurement system is obtained. When the included angle is in the first angle range, the slope inclination of the pixel is determined to be opposite to the sensor observation direction, and the main direction is set to be perpendicular to the slope direction and pointing to the downhill side. When the included angle is in the second angle range, the slope inclination of the pixel is determined to be opposite to the sensor observation direction, and the main direction is set to be perpendicular to the slope direction and pointing to the uphill side. The straight line that passes through the pixel and extends along the main direction is used as the stretching reference line of the pixel.

[0009] As a further aspect of the present invention: the formation process of the contour node sequence is as follows: Select any boundary cell belonging to the target feature from the compensation grid as the starting tracking point. Using the starting tracking point as a reference, explore the cells belonging to the target feature in its neighborhood in a clockwise direction. Connect the adjacent target feature cells in sequence to form a tracking path until return to the starting tracking point and close to form the initial boundary loop. At each connection point of the tracking path, query the stretching amplitude of the pixels on both sides of the connection point in the compensation data set, and calculate the local difference value of the stretching amplitude on both sides. When the local difference value is greater than the preset difference threshold, the location of the connection is marked as a key contour point, and the larger stretching amplitude of the pixels on both sides of the key contour point is assigned to the key contour point as its stretching attribute. All key contour points are arranged in the order of the tracking path to form a sequence of contour nodes with stretching attributes.

[0010] As a further aspect of the present invention: the calculation of the local difference value of the tensile amplitude on both sides specifically includes: The stretching amplitudes of the pixels on both sides of the connection are respectively recorded as the first amplitude and the second amplitude. Taking the connection as the reference, a preset number of consecutive pixels are taken forward and backward along the extension direction of the tracking path, and the stretching amplitudes of the selected pixels are respectively formed into the first amplitude sequence and the second amplitude sequence. Calculate the local median amplitude of the first amplitude sequence and the second amplitude sequence respectively, and subtract the second amplitude from the first amplitude to obtain the amplitude deviation; Divide the amplitude deviation by the global maximum amplitude value of all amplitudes in the first and second amplitude sequences to obtain the normalized deviation ratio, and use the normalized deviation ratio as the local difference value at the connection point.

[0011] As a further aspect of the present invention: the dynamic determination of the simplification threshold of each node based on the stretching properties of each node in the contour node sequence specifically includes: Extract the current node, its preceding neighbor, and its following neighbor from the contour node sequence. Then, combine the stretch attribute values ​​of the preceding neighbor, the current node, and the following neighbor into a local attribute triplet according to the node arrangement order. Find the attribute span value between the maximum and minimum values ​​in the local attribute triplet, and determine the attribute fluctuation intensity of the current node based on the attribute span value; The intensity of attribute fluctuation is used as the adjustment weight to perform weighted correction on the preset simplification baseline value. The corrected value is used as the simplification threshold of the current node. The greater the intensity of attribute fluctuation, the higher the simplification threshold, so as to preserve the contour node at the point of abrupt change in stretched attributes.

[0012] As a further aspect of the present invention: determining the attribute fluctuation intensity of the current node based on the attribute span value specifically includes: The attribute span value is obtained by subtracting the maximum and minimum values ​​in the local attribute triples, and the attribute span value is used as the initial fluctuation index of the current node. Extract the attribute span values ​​of all nodes before the current node from the contour node sequence to form a forward span sequence, and extract the attribute span values ​​of all nodes after the current node to form a backward span sequence. Calculate the midpoint value of each of the forward and backward span sequences, and take the smaller of the two as the local reference span. The initial volatility index is compared with the local benchmark span. If the initial volatility index is greater than a set multiple of the local benchmark span, the current node's attribute volatility intensity is assigned a high volatility intensity; otherwise, it is assigned a low volatility intensity.

[0013] The beneficial effects of this invention are: (1) By calculating the terrain slope and aspect pixel by pixel based on digital elevation data, and combining the observation angle in the rational polynomial parameters, non-uniform stretching is performed on the pixels belonging to the target features in the binary raster along the aspect direction. This restores the width of the features that were compressed on the image due to the terrain perspective contraction to the horizontal projection width. The compensated raster has sufficient pixel spacing in the slope direction to characterize the original geometric turning features of the features. This allows the initial contour extracted by subsequent boundary tracking to truly reflect the actual shape of the features under the horizontal projection, avoiding the systematic deviation of the contour distortion caused by perspective contraction in the subsequent simplification process, and improving the geometric accuracy of feature boundary extraction in remote sensing images.

[0014] (2) By synchronously acquiring the stretching amplitude of each contour node in the compensation data group during the boundary tracking process, the contour node carries a stretching attribute that can reflect the degree of local terrain undulation. In the subsequent simplification process, the simplification threshold of each node is dynamically determined based on its own stretching attribute. This allows contour nodes with abrupt stretching amplitude changes (corresponding to areas of drastic terrain changes) to obtain a higher simplification threshold and be retained, while nodes with gradual stretching amplitude changes to obtain a lower simplification threshold and be eliminated. This method makes the contour simplification degree adaptable to the spatial distribution of terrain undulation. The simplified contour maintains sufficient simplicity in flat areas, while retaining necessary geometric details in areas of drastic terrain changes, thus balancing the cartographic simplicity of the vector results with the terrain fit. Attached Figure Description

[0015] The invention will now be further described with reference to the accompanying drawings.

[0016] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is a flowchart of the process for generating a binary raster that records the spatial distribution of target features in this invention; Figure 3 This is a flowchart illustrating the composition of the compensation data group in this invention. Detailed Implementation

[0017] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] Please see Figure 1 As shown, this invention is an automatic geographic information extraction method based on remote sensing imagery, comprising the following steps: S1, acquire remote sensing images of the target area, rational polynomial parameters corresponding to the remote sensing images, and digital elevation data of the target area; S2 performs pixel-by-pixel feature attribution determination on remote sensing images and generates a binary raster that records the spatial distribution of target features. S3. Based on digital elevation data, the terrain slope and aspect are calculated pixel by pixel. Combined with the observation angle in the rational polynomial parameters, the pixels belonging to the target features in the binary grid are non-uniformly stretched along the aspect direction to restore the width of the compressed features to the horizontal projection width, thus obtaining a compensation grid and recording the stretching amplitude of each pixel to form a compensation data set. S4, perform boundary tracking on the compensation grid in the compensation data group, extract the initial contour nodes of the target features, and simultaneously obtain the stretching amplitude of each node in the compensation data group to form a sequence of contour nodes with stretching attributes. S5 dynamically determines the simplification threshold of each node based on its own stretching properties in the contour node sequence, removes nodes with vertical distances less than the simplification threshold, obtains the simplified contour, and converts the simplified contour into a geographic coordinate vector.

[0019] In S1, remote sensing images of the target area, the corresponding rational polynomial parameters of the remote sensing images, and digital elevation data of the target area are acquired, specifically including: Based on the pre-defined geographic coordinates of the target area, the spatial area to be covered is determined. The remote sensing imagery was acquired from the Resource Satellite Application Center, using panchromatic and multispectral fusion images taken by the Gaofen-2 satellite, with a spatial resolution better than 2 meters. The image format is GeoTIFF, and the image file package includes a rational polynomial parameter file. This file records geometric data such as sensor attitude angle, orbital position, camera focal length, and line-of-sight angle at the time of satellite imaging in text form, which is used for subsequent image geolocation and projection transformation.

[0020] Digital elevation data covering the same target area was obtained from a geospatial data cloud platform. The data used was the global digital elevation data version generated by the US Space Shuttle Radar Topographic Mapping Mission, with a spatial resolution of 30 meters, an elevation accuracy better than 16 meters, and an Img format.

[0021] After obtaining the above data, the spatial reference is unified for the remote sensing image and the digital elevation data: with the projection parameters attached to the remote sensing image as a reference, the digital elevation data is resampled to the same projection coordinate system as the remote sensing image, the bilinear interpolation method is used to avoid introducing obvious jagged edges, and the two are accurately registered in space to ensure that the row and column number deviation of the same ground feature points does not exceed 0.5 pixels.

[0022] After registration is completed, the registered remote sensing image, the rational polynomial parameter file attached to the remote sensing image, and the registered digital elevation data are stored in the local storage device as the basic data set for subsequent pixel-by-pixel processing.

[0023] If the target area spans the coverage of a single remote sensing image, then multiple adjacent images covering the target area are selected, and each image is orthorectified according to rational polynomial parameters before being stitched together so that the stitched image completely covers the target area. At the same time, the digital elevation data is cropped to a spatial range that is completely consistent with the target area.

[0024] All the above data are stored with their original numerical accuracy intact, without any preprocessing operations that would affect the spectral or elevation values ​​of ground features, in order to ensure the accuracy of the original data for ground feature attribution determination and terrain parameter calculation in subsequent steps.

[0025] Please see Figure 2 As shown, in S2, pixel-by-pixel feature attribution is performed on the remote sensing image, generating a binary raster recording the spatial distribution of target features. Specifically, this includes: Using the registered remote sensing image from step S1 as input, determine the number of spectral bands and the center wavelength of each band.

[0026] For the target ground features to be extracted, at least two spectral bands are selected based on their reflectance differences in different spectral bands. For example, the green light band and the near-infrared band are used as the basic calculation bands, and the digital values ​​of each pixel in the remote sensing image are read in these two bands.

[0027] For each pixel, the difference between the values ​​in the two selected bands is calculated. This difference reflects the degree of deviation of the ground object at that pixel in the selected spectral range. The larger the difference, the more obvious the difference between the probability that the pixel belongs to the target ground object and the probability that it belongs to other ground objects.

[0028] The difference calculated for each pixel is arranged sequentially according to the row and column positions of the pixel in the remote sensing image to form a first intermediate grid with the same spatial size as the remote sensing image.

[0029] For each pixel in the first intermediate grid, the neighborhood window size is preset to 3 rows by 3 columns. Centered on the current pixel, the values ​​of all 9 pixels in the window are extracted to form a local value sequence.

[0030] The dispersion value of the local numerical sequence is calculated as follows: the arithmetic mean of the nine values ​​in the sequence is obtained, the absolute value of the deviation of each value in the sequence from the arithmetic mean is calculated, the nine absolute deviation values ​​are summed and then divided by 9. The resulting value is taken as the dispersion value at the pixel position.

[0031] This dispersion value is used to characterize the local fluctuation of the difference between the current pixel and its neighboring pixels. A preset dispersion threshold is set in advance, which is the arithmetic mean of the dispersion values ​​of all pixels in the first intermediate grid multiplied by 0.8.

[0032] The dispersion value of the current pixel is compared with the preset dispersion threshold: if the dispersion value of the current pixel is greater than the preset dispersion threshold, it indicates that the value at the pixel location fluctuates greatly, and the difference value may be affected by isolated noise. In this case, the original value of the pixel in the first intermediate grid is replaced with the median value of the 9 pixels in its neighborhood window to reduce isolated noise points; if the dispersion value of the current pixel is less than or equal to the preset dispersion threshold, it indicates that the pixel value is relatively consistent with the surrounding pixel values ​​and belongs to a continuous area. In this case, the original value of the pixel in the first intermediate grid is retained unchanged.

[0033] After traversing all pixels of the first intermediate grid and completing the above correction, the second intermediate grid is obtained.

[0034] For the second intermediate raster, a preset binarization threshold is determined. The preset binarization threshold is determined by statistically analyzing the distribution of all pixel values ​​in the second intermediate raster and taking the median value of all pixel values ​​as the preset binarization threshold.

[0035] The value of each pixel in the second intermediate grid is compared one by one with the preset binarization threshold: pixels with values ​​greater than or equal to the preset binarization threshold are marked as target features and assigned a value of 1; pixels with values ​​less than the preset binarization threshold are marked as background and assigned a value of 0. After all pixels are marked, the binary grid is generated. The spatial size of the binary grid is consistent with the remote sensing image, wherein the positions of pixels with a value of 1 constitute the spatial distribution range of the target features, and the positions of pixels with a value of 0 constitute the background area.

[0036] Please see Figure 3 As shown, in S3, the terrain slope and aspect are calculated pixel by pixel based on digital elevation data. Combined with the observation view angle in the rational polynomial parameters, the pixels belonging to the target feature in the binary raster are non-uniformly stretched along the slope direction to restore the width of the compressed feature to the horizontal projection width. This results in a compensated raster, and the stretching amplitude of each pixel is recorded, forming a compensated data set, specifically including: Read the elevation values ​​of all pixels in the digital elevation data that are consistent with the spatial range of the binary raster.

[0037] For the cell to be processed, the elevation value of the cell in the digital elevation data is obtained and recorded as the first elevation value. The elevation values ​​of the adjacent cells to the east, west, north, and south of the cell are recorded as the second, third, fourth, and fifth elevation values, respectively.

[0038] Calculate the difference between the second and third elevation values ​​to obtain the horizontal elevation difference; calculate the difference between the fourth and fifth elevation values ​​to obtain the vertical elevation difference. Determine the horizontal ground distance corresponding to one pixel of the digital elevation data, denoted as the first ground distance. Divide the horizontal elevation difference by the first ground distance to obtain the horizontal elevation change rate; divide the vertical elevation difference by the first ground distance to obtain the vertical elevation change rate. The expression for calculating the slope angle is: ; in, Indicates the slope angle, This represents the second elevation value. This represents the third elevation value. This represents the fourth elevation value. This represents the fifth elevation value. Indicates the first ground distance. To represent the square root operation, This indicates the arctangent operation.

[0039] Divide the vertical elevation change rate by the horizontal elevation change rate, and take the arctangent of the resulting ratio to obtain the aspect angle, which represents the azimuth angle of the slope at that pixel. Repeat this process to iterate through all pixels in the digital elevation data corresponding to the locations of all pixels belonging to the target feature in the binary raster, obtaining the slope angle and aspect angle at each target feature pixel.

[0040] The observation angle during satellite transit is read from the rational polynomial parameters. This observation angle represents the azimuth angle projected onto the horizontal plane by the sensor's observation direction. Using the current pixel location as a reference, the values ​​of the slope angle and the observation angle are mapped to a circular metric system ranging from 0 to 360 degrees. The angle between the two values ​​in this circular metric system is then calculated, with the angle ranging from 0 to 180 degrees.

[0041] The first angle range is preset to 0 degrees to 90 degrees, and the second angle range is preset to 90 degrees to 180 degrees.

[0042] When the included angle is within the first angle range, it is determined that the slope inclination of the pixel and the sensor observation direction are in opposite positions on the horizontal projection. In this case, the principal direction is set to be perpendicular to the slope direction and pointing towards the downhill side. When the included angle is within the second angle range, it is determined that the slope inclination of the pixel and the sensor observation direction are in opposite positions on the horizontal projection. In this case, the principal direction is set to be perpendicular to the slope direction and pointing towards the uphill side. The straight line passing through the pixel and extending along the principal direction is used as the stretching reference line of the pixel.

[0043] Based on the slope angle, the number of grid positions that need to be expanded for the pixel is determined in the vertical direction of the stretching baseline.

[0044] The ground distance corresponding to one pixel of the remote sensing image is determined and denoted as the second ground distance. The number of raster locations to be augmented is calculated using the following expression: ; in, Indicates the number of raster positions that need to be expanded. Indicates the first ground distance. Indicates the slope angle, This represents the numerical value of the slope angle. This represents the numerical value of the observed line of sight angle. Represents the cosine operation. This represents the tangent operation. Indicates the second ground distance. This indicates the floor function. When the calculated result... When the value is less than 1, the pixel's The value is set to 1, meaning that each pixel belonging to the target object retains at least one pixel location.

[0045] According to the number of grid positions The binary attribute value of the current cell is migrated from its original position to the newly added grid position along the stretching baseline. The migration direction is along the stretching baseline pointing downhill, and the migration interval is one cell width. That is, for each cell position migrated, a new grid position is added in that direction, and the attribute value of the newly added position is assigned the same attribute value as the current cell.

[0046] The number of raster positions actually migrated by the current pixel is saved as the stretching amplitude of that pixel to temporary storage space. The stretching amplitude is the calculated value. The process iterates through all pixels belonging to the target feature in the binary raster, sequentially performing the aforementioned slope angle calculation, aspect angle calculation, angle determination, stretching baseline setting, amplification quantity calculation, and attribute value migration operations. After the traversal, all target feature pixels after the migration operation and their newly added attribute values ​​together constitute a compensation raster, the spatial range of which is consistent with the binary raster. Simultaneously, the position of each pixel belonging to the target feature and its corresponding stretching amplitude are saved as a key-value correspondence. The correspondence between all pixel positions and stretching amplitudes forms a compensation data set for subsequent steps.

[0047] In S4, boundary tracing is performed on the compensation rasters in the compensation data set to extract the initial contour nodes of the target features. Simultaneously, the stretching magnitude of each node in the compensation data set is acquired, forming a sequence of contour nodes with stretching attributes. Specifically, this includes: In the compensation grid, each pixel is scanned row by row in the order from top to bottom and from left to right, and the first pixel that belongs to the target feature and whose four neighboring areas contain at least one background feature is determined as the starting tracking point.

[0048] Using the starting tracking point as a reference, and taking it as the current tracking point, boundary tracking is performed using an eight-neighborhood clockwise tracking rule: starting from the neighborhood position directly above the current tracking point, the eight neighboring pixels around the current tracking point are explored sequentially in a clockwise direction. The first neighboring pixel belonging to the target feature is identified as the next tracking point, and the connection between the current tracking point and the next tracking point is recorded as a connection point of the tracking path. The next tracking point is updated as the new current tracking point, and the exploration continues in a clockwise direction until the tracking path returns to the starting tracking point and the exploration direction is consistent with the initial direction. When the tracking stops, all connections are closed in the tracking sequence to form the initial boundary loop.

[0049] If, during the tracking process, there are no pixels belonging to the target feature in any of the eight neighborhoods surrounding a current tracking point, then the starting tracking point is abandoned, and the next boundary pixel belonging to the target feature that has not yet been tracked is selected in the compensation grid as the new starting tracking point. Tracking is restarted until all connected regions belonging to the target feature in the compensation grid are traversed, and each connected region forms an initial boundary loop.

[0050] At each connection point of the tracking path, the positions of the two pixels connected at that connection point in the compensation grid are determined and denoted as the first neighboring pixel and the second neighboring pixel, respectively. Using that connection point as a reference, the stretching amplitude corresponding to the first neighboring pixel in the compensation data set is queried and denoted as the first original amplitude; the stretching amplitude corresponding to the second neighboring pixel in the compensation data set is queried and denoted as the second original amplitude.

[0051] Along the extension direction of the tracking path, starting from the connection point, a preset number of consecutive pixels are traced forward, where the preset number is 5 pixels. The stretching amplitudes of the 5 traced pixels in the compensation data group are arranged in the tracing order to form a first amplitude sequence. Along the extension direction of the tracking path, starting from the connection point, a preset number of consecutive pixels are traced backward, where the preset number is also 5 pixels. The stretching amplitudes of the 5 traced pixels in the compensation data group are arranged in the tracing order to form a second amplitude sequence.

[0052] Calculate the median value of the five values ​​in the first amplitude sequence and use this median value as the first local median amplitude; calculate the median value of the five values ​​in the second amplitude sequence and use this median value as the second local median amplitude.

[0053] The difference between the first local median amplitude and the second local median amplitude is taken as the amplitude deviation at that connection. All connections along the tracking path are traversed, and the amplitude deviation at each connection is calculated. The maximum amplitude deviation is then determined as the global maximum amplitude.

[0054] Divide the amplitude deviation at any connection point by the global maximum amplitude value, and the resulting ratio is taken as the local difference value at that connection point. The local difference value ranges from 0 to 1. Calculate the local difference value for each connection point in the tracking path in the same manner.

[0055] A preset difference threshold is set in advance. The preset difference threshold is determined by: statistically analyzing the local difference values ​​at all connections in the tracking path, and taking the arithmetic mean of all local difference values ​​as the preset difference threshold.

[0056] When the local difference value at a certain connection point is greater than the preset difference threshold, the location of the connection point is marked as a key contour point, and the stretching amplitudes of the pixels on both sides of the connection point in the compensation data group are queried. The larger of the two values ​​is assigned to the key contour point as its stretching attribute. When the local difference value at a certain connection point is less than or equal to the preset difference threshold, the connection point is not designated as a key contour point and is not assigned a stretching attribute.

[0057] All key contour points are arranged sequentially according to their order on the tracking path to form a sequence of contour nodes with stretching properties. Each node in the contour node sequence contains the spatial coordinates of the node in the compensation grid and the corresponding stretching property value, which can be called in subsequent steps.

[0058] In S5, the simplification threshold of each node is dynamically determined based on its own stretching properties in the contour node sequence. Nodes with vertical distances less than the simplification threshold are removed to obtain a simplified contour. The simplified contour is then converted into a geographic coordinate vector, specifically including: Read the sequence of contour nodes with stretching properties, where each node contains spatial coordinates and a corresponding stretching property value.

[0059] According to the order of the nodes in the sequence, starting from the second node of the sequence, process them sequentially up to the second to last node.

[0060] For the currently processed node, obtain the stretch attribute value of the node located at the previous position in the contour node sequence as the first attribute value; obtain the stretch attribute value of the node itself as the second attribute value; obtain the stretch attribute value of the node located at the next position as the third attribute value.

[0061] The first attribute value, the second attribute value, and the third attribute value are arranged into local attribute triples according to the node arrangement order.

[0062] Based on the local attribute triplet, the maximum and minimum values ​​among the three values ​​are taken, and the difference between the maximum and minimum values ​​is calculated. This difference is used as the attribute span value of the current node, and the attribute span value is used as the initial fluctuation index of the current node.

[0063] Extract the attribute span values ​​of all nodes preceding the current node from the contour node sequence, and construct a forward span sequence according to the order of these nodes in the sequence; extract the attribute span values ​​of all nodes following the current node from the contour node sequence, and construct a backward span sequence according to the order of these nodes in the sequence.

[0064] Sort all values ​​in the forward span sequence in ascending order, and take the value in the middle position after sorting. If the number of nodes in the forward span sequence is odd, take the value in the middle position as the median value of the forward sequence. If the number of nodes in the forward span sequence is even, take the arithmetic mean of the two middle values ​​as the median value of the forward sequence.

[0065] The same method is used to calculate the median of the backward span sequence. The median of the forward sequence is compared with the median of the backward sequence, and the smaller of the two values ​​is taken as the local reference span of the current node.

[0066] Divide the initial volatility index of the current node by the local benchmark span to obtain the volatility ratio; a multiplier threshold is preset, with a value of 2.0. When the volatility ratio is greater than 2.0, the attribute volatility intensity of the current node is assigned a high volatility intensity; when the volatility ratio is less than or equal to 2.0, the attribute volatility intensity of the current node is assigned a low volatility intensity.

[0067] By traversing all nodes in the contour node sequence except for the first and last nodes in the above manner, the attribute fluctuation intensity of each node is obtained.

[0068] A preset simplified reference value is set in advance. The preset simplified reference value is obtained by taking the ground distance corresponding to one pixel of the compensation grid as the reference distance value, and multiplying the reference distance value by 0.5 as the preset simplified reference value.

[0069] For the current node in the contour node sequence, the preset simplification benchmark value is weighted and corrected based on its attribute fluctuation intensity. Specifically, when the attribute fluctuation intensity of the current node is high, the weighting weight is set to 1.5, and the preset simplification benchmark value is multiplied by 1.5; the product is used as the simplification threshold of the current node. When the attribute fluctuation intensity of the current node is low, the weighting weight is set to 0.8, and the preset simplification benchmark value is multiplied by 0.8; the product is used as the simplification threshold of the current node. The simplification thresholds of all nodes in the contour node sequence, excluding the first and last nodes, are determined sequentially using the above method.

[0070] For each node in the contour node sequence except for the first and last nodes, calculate the vertical distance from the node to the straight line segment connecting its preceding and following adjacent nodes, and use this vertical distance as the vertical distance of the node.

[0071] The vertical distance of the node is compared with the simplification threshold of the node: when the vertical distance is less than the simplification threshold, the node is removed from the contour node sequence; when the vertical distance is greater than or equal to the simplification threshold, the node is retained in the contour node sequence.

[0072] After processing all nodes except the first and last nodes in sequence, the remaining nodes form a simplified contour node sequence in their original order. The first and last nodes in this simplified contour node sequence are the same as the first and last nodes in the contour node sequence.

[0073] The spatial coordinates of each node in the simplified contour node sequence are geocoded and projected according to the rational polynomial parameters described in step S1 to obtain the latitude and longitude coordinates or projected coordinates of each node in the target geographic coordinate system. The transformed coordinates are then connected sequentially according to the node order to form a simplified contour result in the form of geographic coordinate vectors, which can be used for subsequent spatial analysis and mapping applications.

[0074] The working principle of this invention is as follows: Remote sensing images of the target area and their accompanying rational polynomial parameters are acquired, along with digital elevation data for the same area. These are then registered to a unified spatial reference. A pixel-by-pixel feature attribution determination is performed on the registered remote sensing images, generating a binary raster recording the spatial distribution of the target features. Based on the digital elevation data, the terrain slope and aspect are calculated pixel by pixel. Combining this with the observation angle in the rational polynomial parameters, non-uniform stretching is applied to the pixels belonging to the target features in the binary raster along the slope direction. The stretching amplitude of each pixel is determined based on the spatial angular relationship between the slope magnitude, slope aspect, and observation angle, ensuring that the affected area is uniformly stretched. The width of features compressed due to terrain perspective contraction is restored to the horizontal projection width, resulting in a compensation raster and a compensation data set containing the stretching amplitude of each pixel. Boundary tracing is then performed on the compensation raster to extract the initial contour nodes of the target features, and the stretching amplitude corresponding to each node is obtained by querying the compensation data set simultaneously, forming a sequence of contour nodes with stretching attributes. Finally, the simplification threshold of each node is dynamically determined based on its own stretching attribute in the contour node sequence, and nodes with vertical distances less than the simplification threshold are removed to obtain simplified contours. The simplified contours are then converted into geographic coordinate vectors to complete automatic geographic information extraction.

[0075] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. An automatic geographic information extraction method based on remote sensing imagery, characterized in that, Includes the following steps: S1, acquire remote sensing images of the target area, rational polynomial parameters corresponding to the remote sensing images, and digital elevation data of the target area; S2 performs pixel-by-pixel feature attribution determination on remote sensing images and generates a binary raster that records the spatial distribution of target features. S3. Based on digital elevation data, the terrain slope and aspect are calculated pixel by pixel. Combined with the observation angle in the rational polynomial parameters, the pixels belonging to the target features in the binary grid are non-uniformly stretched along the aspect direction to restore the width of the compressed features to the horizontal projection width, thus obtaining a compensation grid and recording the stretching amplitude of each pixel to form a compensation data set. S4, perform boundary tracing on the compensation rasters in the compensation data group, extract the initial contour nodes of the target features, and simultaneously obtain the stretching amplitude corresponding to each node in the compensation data group to form a contour node sequence with stretching attributes. The formation process of the contour node sequence is as follows: Select any boundary cell belonging to the target feature from the compensation grid as the starting tracking point. Using the starting tracking point as a reference, explore the cells belonging to the target feature in its neighborhood in a clockwise direction. Connect the adjacent target feature cells in sequence to form a tracking path until return to the starting tracking point and close to form the initial boundary loop. At each connection point of the tracking path, query the stretching amplitude of the pixels on both sides of the connection point in the compensation data set, and calculate the local difference value of the stretching amplitude on both sides. When the local difference value is greater than the preset difference threshold, the location of the connection is marked as a key contour point, and the larger stretching amplitude of the pixels on both sides of the key contour point is assigned to the key contour point as its stretching attribute. All key contour points are arranged in the order of the tracking path to form a sequence of contour nodes with stretching attributes. The calculation of the local difference in the tensile amplitude on both sides specifically includes: The stretching amplitudes of the pixels on both sides of the connection are respectively recorded as the first amplitude and the second amplitude. Taking the connection as the reference, a preset number of consecutive pixels are taken forward and backward along the extension direction of the tracking path, and the stretching amplitudes of the selected pixels are respectively formed into the first amplitude sequence and the second amplitude sequence. Calculate the local median amplitude of the first amplitude sequence and the second amplitude sequence respectively, and subtract the second amplitude from the first amplitude to obtain the amplitude deviation; Divide the amplitude deviation by the global maximum amplitude value of all amplitudes in the first and second amplitude sequences to obtain the normalized deviation ratio, and use the normalized deviation ratio as the local difference value at the connection point. S5 dynamically determines the simplification threshold of each node based on its own stretching properties in the contour node sequence, removes nodes with vertical distances less than the simplification threshold, obtains the simplified contour, and converts the simplified contour into a geographic coordinate vector.

2. The automatic geographic information extraction method based on remote sensing imagery according to claim 1, characterized in that, The generation of the binary raster recording the spatial distribution of target features specifically includes: Select the values ​​of the target ground features in at least two spectral bands in the remote sensing image, calculate the value difference of the same pixel in different spectral bands, and arrange the obtained differences according to the pixel position to form the first intermediate grid indicating the degree of difference between the target ground features and other ground features in spectral response. For each pixel in the first intermediate grid, the values ​​of all pixels in its neighborhood are taken to form a local numerical sequence. The dispersion value of the local numerical sequence is calculated and compared with a preset dispersion threshold. Based on the comparison result, the value of the corresponding pixel in the first intermediate grid is corrected to weaken isolated noise pixels and retain continuous edge regions, thus forming the second intermediate grid. The value of each pixel in the second intermediate grid is compared with a preset binarization threshold. Pixels that are greater than or equal to the preset binarization threshold are marked as target features, and pixels that are less than the preset binarization threshold are marked as background, thus generating a binary grid.

Images