Terrain Recognition and 3D Reconstruction Method Based on Remote Sensing Mapping Images

Abstract

The present invention relates to the field of remote sensing mapping technology, specifically a method for terrain recognition and three-dimensional reconstruction based on remote sensing mapping images, comprising the following steps: based on remote sensing images, identify and select key spectral bands, adjust the spectral filter configuration, extract key terrain data according to the configuration data, utilize the terrain feature data, analyze terrain changes, record relevant bands and positions, form a comparison table, and generate a three-dimensional model of the terrain by combining natural color and surface texture analysis. In the present invention, by using the optimized band index to screen the key spectral sections, the specificity and efficiency of data acquisition are improved. By precisely adjusting the spectral response interval, it is ensured that key attention is paid to the spectral regions with significant surface reflection characteristics, thereby improving the practical value of spectral data. The spectral configuration data optimizes the response performance of the sensor, making the data collected from the terrain more sensitive and accurate.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing mapping, and particularly to a terrain recognition and three-dimensional reconstruction method based on remote sensing mapping images. Background Art

[0002] Remote sensing mapping technology is a technology that uses remote sensing equipment to collect information on the Earth's surface from a distance, and is widely used in fields such as geology, agriculture, forestry, hydrology, environmental monitoring, and urban planning. It includes using aircraft, satellites, or ground-based sensors to obtain images and other data. After the data is processed, maps, models, or other visual products can be generated to help better understand and manage natural resources.

[0003] Among them, the terrain recognition and three-dimensional reconstruction method of remote sensing mapping images involves using remote sensing images for automatic terrain recognition and three-dimensional model construction. The purpose is to identify specific terrain features, such as mountains, rivers, buildings, etc., through high-resolution remote sensing images, and convert two-dimensional data into three-dimensional models. The three-dimensional models are of great importance for performing geospatial analysis, creating disaster management predictions, etc., and can greatly improve the visualization and analysis capabilities of terrain data.

[0004] Existing technologies often fail to accurately distinguish terrain details due to insufficient use of spectral data and limitations of processing algorithms when dealing with areas with complex or rapidly changing terrains, such as rugged mountainous areas or dense urban environments. The generalization processing of spectral data and the reliance on simplified stereo matching techniques often perform poorly in application scenarios with high-precision processing requirements, and can easily lead to serious consequences in urban planning, disaster management, or geospatial analysis that require highly accurate data support, such as planning mistakes or insufficient early warning responses, affecting decision-making quality and safety management efficiency. Summary of the Invention

[0005] The purpose of the present invention is to solve the disadvantages existing in the prior art, and to propose a terrain recognition and three-dimensional reconstruction method based on remote sensing mapping images.

[0006] To achieve the above purpose, the present invention adopts the following technical solution: A terrain recognition and three-dimensional reconstruction method based on remote sensing mapping images, including the following steps:

[0007] S1: Based on the surface reflection spectral data in the remote sensing mapping image, record the spectral intensity under different bands, identify the bands with abnormal spectral intensity according to the band number and light intensity, and summarize them as the preferred band index;

[0008] S2: Based on the preferred band index, adjust the spectral filter configuration, read the current configuration parameters of the hyperspectral sensor, test the adjustable range of the spectral response section, match the main response section of each band, retain the intersection section with the preferred band index, and reset the cut-off band and central band of the spectral filter to generate spectral configuration data;

[0009] S3: According to the spectral configuration data, extract key terrain data from the slope information layer, curvature raster layer, and elevation gradient layer, calculate the ratio change of slope and curvature, identify the spatially continuous section, then number the section and map its geometric features to obtain the terrain feature geometry map;

[0010] S4: According to the abnormal change area in the terrain feature geometry map, analyze the bands with weak response associated with the spectral configuration data, record the central position of the area, and pair each position with the band number, and integrate them into a position-band comparison table.

[0011] The improvements of the present invention are that the preferred band index includes spectral intensity anomaly index, spectral adjustment reaction, and band feature identifier, the spectral configuration data includes band adjustment range, filter characteristics, and spectral response configuration, the terrain feature geometry map includes terrain continuity identification result, geometric feature mapping, and spatial change analysis result, and the position-band comparison table includes central position mark, band response matching result, and terrain anomaly index.

[0012] The improvements of the present invention are that the specific steps for obtaining the preferred band index are as follows:

[0013] S111: Based on the surface reflection spectral data in the remote sensing mapping image, identify the difference in illumination intensity and spectral intensity of each band according to the band number and spectral intensity, construct the mapping relationship between the band number and spectral intensity, and form a spectral data mapping structure;

[0014] S112: Using the spectral data mapping structure, calculate the absolute difference between the spectral intensity and illumination intensity of each band, identify the bands exceeding the fluctuation threshold, and determine the band number interval with anomalies to obtain the abnormal band interval;

[0015] S113: Through the abnormal band interval, analyze the band numbers within the interval, compare the difference between the spectral intensity and illumination intensity, classify the band numbers according to the band number frequency and difference degree, and summarize them into the preferred band index.

[0016] The improvements of the present invention are that the specific steps for obtaining the spectral configuration data are as follows:

[0017] S211: Based on the preferred band index, adjust the spectral filter configuration to align the current central band, cut-off band, and response section boundary parameters of the spectral filter. By comparing the center and width, screen the band segments that match the preferred band to obtain the band matching section;

[0018] S212: Invoke the band matching section to refine the response boundary range of each adjustable band, using the formula:

[0019] ;

[0020] Obtain the band overlap amount , sort the bands, screen the band segments with high overlap degree to obtain the band overlap degree list, where represents the overlap width between the th band and the th matching section, is the alignment factor between the th band and the th matching section, represents the bandwidth interaction characteristic between the th band and the th matching section, represents the compatibility parameter between the th band and the th matching section, represents the number of matching sections, represents the total number of bands;

[0021] S213: According to the band overlap degree list, detect the central band and cut-off band of each intersection section, retain the intersection section with the preferred band index, and reset the central band and cut-off band of the spectral filter to generate spectral configuration data.

[0022] An improvement of the present invention is that the step of obtaining the geometric map of the terrain feature is specifically as follows:

[0023] S311: According to the spectral configuration data, call the remote sensing mapping image, extract the terrain slope from the slope information layer, match it with the surface curvature extracted from the curvature raster layer, calculate the ratio of the terrain slope to the surface curvature, screen the spatial points that meet the conditions to obtain the slope-curvature ratio screening set;

[0024] S312: According to the slope-curvature ratio screening set, combine the elevation gradient data, analyze the elevation change of each point, compare the consistency between the slope-curvature ratio and the elevation change, identify and number the continuous spatial sections to obtain the spatial continuity section number;

[0025] S313: Based on the spatial continuity section number, use the formula:

[0026] ;

[0027] Calculate the geometric expression value of the spatial form within the calculation section , and map its geometric features to obtain a geometric map of terrain features, where represents the slope change at the th point in the th section, represents the elevation fluctuation at the th point in the th section, represents the curvature trend at the th point in the th section, represents the boundary area of the th section, represents the edge connectivity of the th section, represents the average value of the internal distribution density of the th section, indicates the number of points contained in the th section.

[0028] The improvement of the present invention is that the acquisition steps of the position-band comparison table are specifically as follows:

[0029] S411: According to the abnormal change areas in the geometric map of terrain features, select boundary undulation sections, slope mutation sections and texture fracture sections, collect image coordinate data and spatial position information of the corresponding sections, call spectral configuration data, match image coordinates with band reflectance sequences, and obtain spatial reflectance corresponding data;

[0030] S412: Based on the spatial reflectance corresponding data, analyze the band reflectance sequence, calculate the absolute difference between the band reflectance and the regional average reflectance at each group of coordinate points, screen the band numbers with small reflectance differences, and record the corresponding image coordinate points to obtain the positions of the bands with weak responses;

[0031] S413: Based on the positions of the bands with weak responses, construct a pairing relationship between coordinate points and band numbers, and integrate them into a position-band comparison table.

[0032] The improvement of the present invention also includes:

[0033] S5: Based on the position-band comparison table, extract the pixel set of the associated area, analyze the natural color layer and surface texture layer of the area, compare the texture color samples by calculating the RGB color mean and local fluctuations, and simulate the terrain surface features to perform three-dimensional reconstruction of the terrain to obtain a three-dimensional terrain model;

[0034] The three-dimensional terrain model includes three-dimensional coordinate mapping, terrain surface simulation information, and hierarchical structure display.

[0035] The improvement of the present invention is that the acquisition steps of the three-dimensional terrain model are specifically as follows:

[0036] Based on the position-band look-up table, extract the pixel set of the associated area, read the RGB channel data of the natural color layer and the surface texture layer, calculate the average value of the RGB channels respectively, and obtain the regional channel color reference;

[0037] Based on the regional channel color reference, analyze the color change between adjacent pixels in each area, and use the formula:

[0038] ;

[0039] Calculate the color volatility , simulate the terrain surface features, perform three-dimensional terrain reconstruction, and obtain the three-dimensional terrain model, where is the number of pixels in the area, , and are the red, green, and blue channel data of the th pixel in the area, , and are the average color data of the RGB channels of the area.

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

[0041] In the present invention, by using the preferred band index to screen the key spectral sections, the specificity and efficiency of data acquisition are improved. By precisely adjusting the spectral response range, it is ensured that key attention is paid to the spectral regions with significant surface reflection characteristics, improving the practical value of spectral data. The spectral configuration data optimizes the response performance of the sensor, making the data collected from the terrain more sensitive and accurate. The detailed construction of the terrain feature geometric map enables the construction of the three-dimensional model to no longer rely only on the simple visual data of the surface layer, but to penetrate into the microscopic structure of the terrain, providing a more dynamic and realistic three-dimensional view. Especially in geographical spatial analysis and disaster prediction, the effectiveness of decision-making is significantly improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 is the main step flow chart proposed by the present invention;

[0043] Figure 2 is the flow chart for obtaining the preferred band index in the present invention;

[0044] Figure 3 This is the flowchart for obtaining the spectral configuration data in the present invention;

[0045] Figure 4 This is the flowchart for obtaining the topographic feature geometry map in the present invention;

[0046] Figure 5 This is the flowchart for obtaining the position - band comparison table in the present invention;

[0047] Figure 6 This is the flowchart for obtaining the three - dimensional topographic model in the present invention;

[0048] Figure 7 This is the topographic feature map of the present invention;

[0049] Figure 8 This is the topographic geometric feature map of the present invention;

[0050] Figure 9 This is the topographic model map of the present invention. Detailed implementation manners

[0051] In order to make the objectives, technical solutions and advantages of the present invention more clear and understandable, the following further details the present invention in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0052] In the description of the present invention, it should be understood that the orientation or positional relationships indicated by the terms "length", "width", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc. are based on the orientation or positional relationships shown in the accompanying drawings. These are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation. Therefore, it should not be construed as a limitation to the present invention. In addition, in the description of the present invention, "a plurality of" means two or more unless otherwise specifically defined.

[0053] Embodiment:

[0054] Please refer to Figure 1 , the present invention provides a technical solution: a method for terrain recognition and three - dimensional reconstruction based on remote sensing mapping images, including the following steps:

[0055] S1: Based on the surface reflection spectral data in the remote sensing mapping image, record the spectral intensities under different bands (including standard bands such as blue light, green light, red light, near - infrared light, etc.). For each band, identify the bands with abnormal spectral intensities according to the band number and light intensity, and summarize them as the preferred band indexes;

[0056] S2: Based on the preferred band index, adjust the spectral filter configuration, read the current configuration parameters of the hyperspectral sensor (including central wavelength, bandwidth, gain, etc.), test the adjustable range of the spectral response section, match the main response section of each band, retain the intersection section with the preferred band index, and reset the cut-off band and central band of the spectral filter to generate spectral configuration data;

[0057] S3: According to the spectral configuration data, call the remote sensing mapping image, extract key terrain data from the slope information layer, curvature raster layer and elevation gradient layer, including terrain slope, curvature and elevation change, calculate the ratio change of slope and curvature, identify the spatially continuous section, then number the section and map its geometric features to obtain the terrain feature geometry map;

[0058] S4: Analyze the response-weakening bands (such as the near-infrared band of water bodies or the short-wave infrared band of vegetation) associated with the spectral configuration data according to the abnormal change areas in the terrain feature geometry map, record the central positions of the areas (using the WGS84 coordinate system), and pair each position with the band number, and integrate them into a position-band comparison table;

[0059] S5: Based on the position-band comparison table, extract the pixel set of the associated area, analyze the natural color layer (RGB band combination) and surface texture layer (calculated by the gray-level co-occurrence matrix GLCM) of the area, calculate the RGB color mean and local fluctuation by calculating, compare the texture color samples, and simulate the terrain surface features to perform three-dimensional reconstruction of the terrain to obtain a three-dimensional terrain model.

[0060] The preferred band index includes spectral intensity anomaly indicators, spectral adjustment responses, and band feature identifiers. The spectral configuration data includes band adjustment ranges, filter characteristics, and spectral response configurations. The terrain feature geometry map includes terrain continuity identification results, geometric feature mappings, and spatial change analysis results. The position-band comparison table includes central position markers, band response matching results, and terrain anomaly indices. The three-dimensional terrain model includes three-dimensional coordinate mappings, terrain surface simulation information, and hierarchical structure displays.

[0061] Spectral intensity refers to the intensity of the electromagnetic wave energy reflected back by different materials on the earth's surface in remote sensing mapping images. The energy exhibits different intensities in different bands (such as visible light, infrared band, etc.), and is usually used to distinguish the materials, components, etc. of surface objects; the optimal band index analyzes the spectral intensities of different bands, selects those bands that can significantly exhibit the characteristics of surface objects, and summarizes the band numbers to form the optimal band index. This process needs to be adjusted according to the lighting conditions and surface reflection characteristics to ensure that the selected bands can effectively identify the characteristics of different surface objects; the spectral filter is configured in the hyperspectral sensor. The spectral filter is used to select and control which bands can be received by the sensor. Adjusting the spectral filter configuration means resetting its working band range and response interval to adapt to specific remote sensing data requirements; the terrain slope, curvature, and elevation change are the basic terrain features describing the surface undulation. The slope represents the inclination degree of the surface, the curvature describes the degree of surface bending, and the elevation change reflects the height difference of the ground undulation, which are the basic data for terrain recognition and analysis; the location-band comparison table is a table that maps specific geographical locations to their corresponding band numbers, and is used to analyze the relationship between the spectral characteristics and terrain characteristics of a specific area; as shown in Table 1, the location-band comparison table is a data that pairs the position coordinates of the spatially abnormal areas identified in the remote sensing image with the band numbers showing weak response in the area. The table contains: spatial coordinates (WGS84, the central geographical location coordinates (latitude and longitude) of each abnormal terrain point or area), image coordinate points (the pixel numbers or row-column positions in the original remote sensing image), the band numbers with weak response (the bands with the weakest spectral response or the smallest difference (such as near-infrared, short-wave infrared)), and the band reflectance values (the reflectance of each weak response band at this point).

[0062] Table 1: Location-Band Comparison Table:

[0063]

[0064] Please refer to Figure 2 , and the specific steps for obtaining the optimal band index are as follows:

[0065] S111: Based on the surface reflection spectral data in the remote sensing mapping image, according to the band number and spectral intensity, identify the difference between the lighting intensity and spectral intensity of each band, construct the mapping relationship between the band number and spectral intensity, and form the spectral data mapping structure;

[0066] Based on the surface reflectance spectral data in the remote sensing mapping images, the reflectance of each pixel at different bands is extracted. The band numbers obtained are set as B1 to B7, and the corresponding reflectances are 0.12, 0.15, 0.13, 0.24, 0.41, 0.36, and 0.39 in sequence. The light intensities are 0.14, 0.13, 0.15, 0.22, 0.43, 0.33, and 0.35 respectively. On this basis, a mapping relationship between the number and the spectral intensity is established. The reflectance data and the light intensity are combined. For each band number, the difference between its spectral intensity and the light intensity is calculated respectively, and the absolute difference method is used for processing, that is, D1 = |0.12 - 0.14| = 0.02, D2 = |0.15 - 0.13| = 0.02, D3 = 0.02, D4 = 0.02, D5 = 0.02, D6 = |0.36 - 0.33| = 0.03, D7 = |0.39 - 0.35| = 0.04. To set the fluctuation benchmark, it is necessary to analyze based on historical sampling data. Assume that the above-mentioned difference data set is obtained from 200 sample points in the early stage, its mean value is 0.025, and the standard deviation is 0.012. Therefore, the fluctuation threshold can be set to 0.025 + 0.012 = 0.037, which is used as the baseline for judging the single-point band difference. Substituting the above D1 to D7 for comparison, it is found that none of them exceed the threshold, only D7 is close to the threshold. The data structure is sorted into a mapping structure recorded by band number, and items such as number, spectral intensity, light intensity, and difference parameters are retained in the content. A spectral data mapping structure is constructed as the input basis for subsequent judgment and screening. This structure can be used for single-point inversion analysis or area comparison analysis of any pixel, providing complete data information support for subsequent fluctuation anomaly screening.

[0067] S112: Using the spectral data mapping structure, calculate the absolute difference between the spectral intensity and the light intensity of each band, identify the bands that exceed the fluctuation threshold, and determine the band number interval where anomalies exist to obtain the abnormal band interval;

[0068] Perform differential threshold judgment for each band number. Set the reference fluctuation threshold to 0.037. Call the aforementioned differential data D1 to D7 one by one, directly compare the differential value of each band number with the threshold, and determine whether it exceeds the reference. After the comparison, it is found that D7 = 0.04 exceeds the threshold, and the others do not exceed the limit. Therefore, B7 is initially identified as an abnormal band. Perform a batch judgment operation on the entire image. Based on all pixels, accumulate the number of times each band number exceeds the limit in each pixel. Set the standard that a certain band exceeding 80% in all pixels is a significant abnormality standard. Assume that among a total of 200 pixels, B7 exceeds the limit in 182 pixels, accounting for 91%, which is much higher than the set standard. Confirm it as an abnormal band. At the same time, perform the same type of judgment on the adjacent band B6. Its number of times exceeding the limit is 120, accounting for 60%, which does not constitute a significant abnormality. However, because B6 is an adjacent band to B7, it needs to be included in the analysis scope and classified into the abnormal band number range. Determine the range of this interval as B6 to B7. This identification process needs to be stably reproduced in the regional distribution, that is, only the bands that maintain the same abnormal pattern in multi-temporal or multi-regional images can be determined as credible abnormal sections. Integrate the number information, differential amplitude, and repeated distribution status to obtain the abnormal band interval.

[0069] S113: Through the abnormal band interval, analyze the band numbers within the interval, compare the difference between the spectral intensity and the illumination intensity, classify the band numbers according to the number frequency and the degree of difference of the bands, and summarize them into the preferred band index;

[0070] Through the abnormal band interval, that is, B6 to B7, reprocess the distribution frequency and the degree of difference of each pixel in the target area. Compare whether the recurrence frequency of this number in the overall sample exceeds 90%. If a certain band number appears stably in all sample areas or multi-temporal data and is accompanied by an excessive difference, it is identified as the target preferred band. When performing such screening, the number frequency and the spectral offset amount need to be used as the core parameters. Call the previous calculation results that the offset of B6 is 0.03 and B7 is 0.04, which respectively represent their actual difference values under the current pixel. At the same time, quote the occurrence ratios of these two numbers in the entire sample area as: B6 is 78% and B7 is 92%. Taking the scoring reference of 0.03 as the standard, exclude B6 and include B7 in the preferred number set. The set scoring model does not use weights but uses the rule that both the difference value and the frequency exceed the reference as the judgment rule. Finally, record B7 as the effective preferred band number and add it to the final index set, and summarize and generate the preferred band index. This index can be called in subsequent application links such as classification recognition, feature recognition, and target extraction.

[0071] Please refer to Figure 3 , the specific steps for obtaining the spectral configuration data are as follows:

[0072] S211: Based on the preferred band index, adjust the spectral filter configuration, align the current central band, cut-off band, and response section boundary parameters of the spectral filter, screen the band segments that match the preferred band by comparing the center and width, and obtain the band matching section;

[0073] Align the current central band, cut-off band, and response section boundary parameters of the spectral filter. First, read the currently set multiple central bands and cut-off bands from the sensor. For example, the current central bands are 520nm, 610nm, 700nm, 820nm, and the corresponding cut-off bands are 500nm, 590nm, 680nm, 800nm. Then, extract the center and cut-off information of the target matching band from the preferred band index. The preferred central bands are set as 525nm, 615nm, 705nm, 815nm, and the cut-off boundaries are set as 505nm, 595nm, 685nm, 805nm. Compare the offset of the center position of each group of bands in turn. For example, the deviation between 520nm and 525nm is 5nm, and the deviation between 610nm and 615nm is 5nm. The deviation values are all controlled within the allowable offset range of ±10nm, and the band pairs with alignable positions can be determined. Subsequently, continue to compare the cut-off bands. For example, 500nm and 505nm, 590nm and 595nm, and the deviations are also within the allowable range, indicating that the cut-off boundaries also have proximity. Then, it is necessary to further examine the response section width of the band, that is, the overall range of each band. For example, for the 520nm central band, its actual spectral response covers 510nm to 530nm. If the corresponding preferred section response range is 515nm to 535nm, then its spectral overlap section can be judged as 515nm to 530nm, the overlap width is 15nm, and the proportion reaches more than 75% of the original bandwidth. This situation can be regarded as having high adaptability. Perform the same operation on other bands in turn. If the overlap section of the 700nm corresponding section is 690nm to 710nm and its coverage rate is 100%, it can be directly judged as a perfect match. Finally, the bands that meet the requirements in terms of center offset, cut-off overlap, and response bandwidth overlap are used as valid matching items, sorted into a unified structure to form the band matching section.

[0074] S212: Call the band matching section, refine the response boundary range of each adjustable band, and use the formula:

[0075] ;

[0076] Obtain the band overlap amount , sort the bands, screen the band segments with high overlap, and obtain the band overlap degree list, where represents the th band and the The overlapping width of the matching segments, where the width refers to the overlapping part of two bands in the spectrum and is used to calculate their overlapping degree. is the alignment factor of the th band and the th matching segment, which quantifies the alignment degree of their central positions. th band and the bandwidth interaction characteristic of the th band and the th matching segment, which is the sum of the total bandwidths of the two bands, representing the overall coverage ability of the two bands, reflecting bandwidth similarity or interference, and measuring whether the spectral coverage overlap of the two bands in the same region is reasonable, used to control whether the two bandwidths are too wide and the matching degree is inflated, and the actual response may be very weak. compatibility parameter of the th band and the th matching segment, which reflects the spectral overlap degree and physical property compatibility.

[0077] Suppose the centers of the currently adjustable 3 bands are 520nm, 610nm, and 700nm respectively, the centers of the matching segments are 525nm, 615nm, and 705nm respectively, the actual bandwidth of each group of bands is uniformly set to 20nm, the response ranges are [510–530], [600–620], [690–710] respectively, and the response range of the matching segments is [515–535], [605–625], [695–715]. According to the intersection of the response ranges, judge the overlapping width. The overlapping part of band 1 and segment 1 is 515–530, and the overlapping width is 15nm. The overlapping of band 2 and segment 2 is 605–620, and the width is 15nm. The overlapping of band 3 and segment 3 is 695–710, which is also 15nm. Then, quantify the center alignment factor, and use the center difference and the reference range conversion. For example, the difference between 520nm of band 1 and 525nm of segment 1 is 5nm. Under the reference width of 30nm, the converted alignment factor is 0.833. The differences between band 2 and band 3 and their corresponding segments are both 5nm, so the alignment factors are both 0.833. Then, evaluate the bandwidth interaction characteristic. The sum of the two bandwidths can be used as the interaction term, which is both 20 + 20 = 40nm. The compatibility parameter is measured by the difference of the normalized response curves. Suppose the maximum deviation conversion between the three groups of spectra is converted into a spectral difference of 10nm, then the compatibility parameter is set to 10. Substitute the above parameters into the following formula to calculate the overlapping amount of each band:

[0078] ;

[0079] The assignments of each parameter are as follows:

[0080] ;

[0081] ;

[0082] , all the same;

[0083] Substitute and calculate the formula items as follows:

[0084] ;

[0085] , from which the coincidence degree of the three groups of bands is 0.4165. The numerical value indicates that the current bands and their corresponding matching sections are consistent in response width, position alignment, and spectral compatibility. The coincidence degree ranking results are equivalent. Therefore, all three are included in the subsequent screening to form a band coincidence degree list.

[0086] S213: According to the band coincidence degree list, detect the central band and cut-off band of each intersection section, retain the intersection sections with the preferred band index, and reset the central band and cut-off band of the spectral filter to generate spectral configuration data;

[0087] First, extract the central and boundary data of the bands with the top-ranked coincidence degrees. For example, if band 1 is 520 nm, the cut-off is 510 nm, and the corresponding preference is 525 nm and 515 nm, then the central offset is 5 nm and the cut-off offset is 5 nm. Determine whether it meets the reconfiguration standard. Assume the central allowable difference is ±10 nm and the cut-off allowable difference is ±10 nm. Then this band is retained. After that, reset the center and cut-off of the retained bands. The central band takes the value of the centroid of adjacent sections. For example, 520 nm and 530 nm are combined into 525 nm, and the cut-off boundary is processed according to the average value of the boundaries. For example, 510 nm and 515 nm are combined into 512.5 nm. Perform this operation on all retained bands in turn to generate a combined central and cut-off parameter in a unified format, forming spectral configuration data for configuring the spectral filter.

[0088] Please refer to Figure 4 , and the specific steps for obtaining the topographic feature geometry map are as follows:

[0089] S311: According to the spectral configuration data, call the remote sensing mapping image, extract the topographic slope from the slope information layer, match it with the surface curvature extracted from the curvature raster layer, calculate the ratio of the topographic slope to the surface curvature, screen the spatial points that meet the conditions, and obtain the slope-curvature ratio screening set;

[0090] Call the remote sensing mapping image data source. This data has a resolution of 1 meter and covers the reflection band. Load the slope information layer and the curvature raster layer in a superposition manner. Unify the pixel numbers of each spatial point in the image with the raster index of the layer, and extract the slope and curvature data of each point. For example, for the point numbered P001, its slope is 19.3 degrees and the corresponding curvature is 0.023 / m². Then the slope-curvature ratio of this point is 19.3 divided by 0.023, which equals 839.13. Process all the overlapping points of the layers in a sequential traversal manner, calculate their respective slope-curvature ratios, and construct a complete ratio dataset. On this basis, calculate the mean and standard deviation of this dataset. For example, the ratio mean is 470 and the standard deviation is 125. Then set the screening threshold as the mean plus one standard deviation, which is 595. Compare the ratio data according to this threshold. For example, if the ratio of the point numbered P002 is 627.1, it is retained; if the ratio of the point numbered P003 is 412.7, it is excluded. Form a point set with all the points whose ratio is greater than the screening threshold. The screened spatial point set is the slope-curvature ratio screening set.

[0091] S312: According to the slope-curvature ratio screening set, combined with the elevation gradient data, analyze the elevation change of each point, compare the consistency between the slope-curvature ratio and the elevation change, identify and number the continuous spatial sections, and obtain the spatial continuity section numbers;

[0092] According to the points contained in the slope-curvature ratio screening set, call the elevation change data that completely coincides with its coordinates in the elevation gradient layer to analyze the elevation fluctuation trend of each point. Use a continuous comparison method to calculate the elevation difference between adjacent points. For example, for the points numbered P101 and P102, their elevations are 256.7 meters and 257.1 meters respectively. Then the elevation difference between them is 0.4 meters. If the judgment condition is set that the elevation difference between continuous points shall not exceed 0.5 meters, then this point pair meets the condition and is used as a candidate continuous point. On this basis, analyze the change trend of the slope-curvature ratio of adjacent points. Use the relative change rate calculation method. If the ratios of two consecutive points are 839.13 and 841.77, then the change rate is the absolute value of (841.77 minus 839.13) divided by 839.13, which is 0.0031, that is, 0.31%. When the change rate is lower than the 5% threshold, it can be judged that the trends are consistent. The judgment standard is that if both the ratio change rate and the elevation difference meet the conditions, they are continuous spatial points. Mark and group them sequentially and number them in an increasing numbering method to construct a section identifier. For example, the numbered section Z01 contains the points P100 to P118, and the section numbers of the continuous point set are output.

[0093] S313: Based on the spatial continuity section numbers, use the formula:

[0094] ;

[0095] Calculate the geometric expression value of the spatial form within the section , and map its geometric features to obtain a geometric map of terrain features, where represents the slope change at the th point in the th section, indicating how the slope of this point changes with the terrain, represents the elevation fluctuation at the th point in the th section, describing the change in the relative height of this point, represents the curvature trend at the th point in the th section, reflecting the curvature change of the terrain at this point, represents the boundary area of the th section, used to describe the total area size of this section, represents the edge connectivity of the th section, describing the spatial connectivity between this section and the surrounding area, represents the mean of the internal distribution density of the th section, describing the distribution density of points within this section and reflecting the uniformity of spatial distribution, represents the number of points contained in the th section;

[0096] Let the section be , which contains a total of points, and the original geometric parameters are as follows:

[0097] Slope change : [22.4, 23.1, 22.9, 21.7, 23.4]°, normalized by the maximum value of 23.4;

[0098] Elevation fluctuation : [15.7, 16.3, 15.9, 15.2, 16.0] m, normalized by the maximum value of 16.3;

[0099] Curvature trend : [0.019, 0.021, 0.020, 0.018, 0.022] 1 / m², normalized by the maximum value of 0.022.

[0100] The corresponding normalized values are:

[0101] ;

[0102] ;

[0103] ;

[0104] Section parameters (original value & normalized value):

[0105] Boundary area , with 1000 as the upper limit of normalization reference;

[0106] Edge connectivity , normalized with 5 as the maximum reference value;

[0107] Average distribution density , normalized with 3 as the maximum reference value.

[0108] Normalization result:

[0109] ;

[0110] ;

[0111] ;

[0112] Substitute into the calculation of the numerator part:

[0113] ;

[0114] ;

[0115] ;

[0116] ;

[0117] ;

[0118] Sum of numerators ;

[0119] Calculate the denominator part:

[0120] ;

[0121] ;

[0122] Denominator ;

[0123] Result:

[0124] ;

[0125] is 7.375, representing the combined strength of the three geometric attributes of slope change, elevation fluctuation, and curvature trend in this section. The larger this value, the more active the geometric characteristics of the spatial points in this section are as a whole, the more obvious the structural differences are, and the relatively higher the comprehensive morphological complexity is. The smaller the value, the flatter the terrain change or the simpler the structure. That is, as a graphic expression factor of a segment, it is mapped as a primitive color or classification number in the terrain feature geometry map. The further arrangement of the expression value can be done by using the natural discontinuity method or the equidistant grading method to classify all the segments. The values ​​are divided into several levels, thus completing the construction of the spatial zoning map of terrain features.

[0126] See also Figure 5 , the steps for obtaining the position band comparison table are as follows:

[0127] S411: According to the abnormal change area in the terrain feature geometry map, the boundary undulation section, the slope mutation section and the texture fracture section are selected, the image coordinate data and spatial position information of the corresponding section are collected, the spectrum configuration data is called, the image coordinates and the band reflectance sequence are matched, and the corresponding data of the spatial reflectance are obtained;

[0128] First, according to the characteristics of boundary undulation segments, slope mutation segments and texture fracture segments, feature segments are selected and their image coordinate data and spatial position information are collected. For each pixel point in the abnormal area, a preliminary analysis of the grayscale value and slope is first performed to mark the area with the most significant change. For example, pixel line segments with a grayscale difference threshold greater than 10 are set as candidates. Then, further classification is performed in combination with the slope mutation or texture fracture characteristics, and the central coordinate points of each type of area are selected for spatial positioning conversion and their coordinate information is recorded. For example, the spatial coordinates of the selected image coordinate point P01 are X=135.2 meters, Y=87.4 meters, and the spatial coordinates of the point P02 are X=142.7 meters, Y=93.2 meters. The information is used to associate the reflectance values ​​of each band in the spectral configuration data, and the band reflectance at each coordinate point is corresponded to its spatial coordinates one by one to form spatial reflectance corresponding data, that is, the spatial coordinates and band reflectance information are combined to establish the corresponding relationship between image coordinates, spatial coordinates and reflectance sequence.

[0129] S412: Based on the spatial reflectivity corresponding data, analyze the band reflectivity sequence, calculate the absolute difference between the band reflectivity at each group of coordinate points and the regional average reflectivity, select the band numbers with small reflectivity differences, and record the corresponding image coordinate points to obtain the weak response band positions;

[0130] Calculate the absolute difference between the band reflectance values of each set of coordinate points and the average reflectance of the region. Take the criterion that the absolute difference is less than a certain threshold (e.g., 0.05) to screen out the band number with the smallest reflectance difference. For example, assume the reflectance of point P01 is B1: 0.42, B2: 0.38, B3: 0.45, and the average reflectance of the region is B1: 0.40, B2: 0.35, B3: 0.43. The calculated absolute differences are B1: 0.02, B2: 0.03, B3: 0.02 respectively. Thus, it is determined that the differences in the B1 and B3 bands are smaller and meet the screening conditions. Finally, record the B1 and B3 bands corresponding to point P01, and correspond this information with the image coordinate points to obtain the position of the band with weak response, that is, screen out the band with the smallest difference in band reflectance and its corresponding coordinates.

[0131] S413: Based on the position of the band with weak response, construct the pairing relationship between the coordinate points and the band numbers, and integrate them into a position-band comparison table;

[0132] Integrate each image coordinate point with its corresponding band number with weak response one by one to form a data set and export it. For example, pair point P01 with the bands B1 and B3 with the smallest reflectance difference to get the corresponding relationship (P01, B1: 0.42, B3: 0.45). Similarly, pair other coordinate points, and integrate all coordinate points and band numbers into a complete position-band comparison table. The table contains each image coordinate point and its corresponding weakest response band, providing basic data for subsequent analysis and processing.

[0133] Please refer to Figure 6 , the steps for obtaining the three-dimensional terrain model are specifically as follows:

[0134] Based on the position-band comparison table, extract the pixel set of the associated region, read the RGB channel data of the natural color layer and the surface texture layer, and calculate the average value of the RGB channels respectively to obtain the regional channel color reference;

[0135] Extract the pixel set of the associated area, which involves the RGB channel data of the natural color layer and the surface texture layer. During the execution process, the data set of each layer is precisely divided, and the color values of the R, G, and B channels are read. The data is extracted and stored through a specific area division mode. Then, statistical processing is performed on the RGB channels of each area. For example, calculate the average color values of the R channel, G channel, and B channel of the pixels in each area. The calculation method is to perform a weighted sum of the color data of all pixels in each area and then divide by the total number of pixels to obtain the color benchmark of the area. For example, for a certain area, if this area contains 1000 pixels and the pixel values of its R channel are [120, 130, 110,...], the calculated average R channel color is 115. Similarly, the average values of the G and B channels will also be calculated based on all the pixels in each area. The obtained average RGB channel values of the area are the color benchmarks of the area.

[0136] Based on the area channel color benchmark, analyze the color changes between adjacent pixels in each area, using the formula:

[0137] ;

[0138] Calculate the color volatility , which represents the overall fluctuation degree of the pixel colors in this area, simulate the terrain surface features, perform 3D terrain reconstruction, and obtain the 3D terrain model. Among them, is the number of pixels in the area, , and are the red, green, and blue channel data of the th pixel in the area, , and are the average color data of the RGB channels of the area;

[0139] First, sum the color difference values of each pixel, and then average by the number of pixels to obtain the color volatility index of this area. If there are 3 pixels in the y area, and the values of the R, G, and B channels are as follows:

[0140] For pixel 1: ;

[0141] For pixel 2: ;

[0142] For pixel 3: ;

[0143] The average values (benchmark values) of the RGB channels of area y are:

[0144] , , ;

[0145] For pixel 1, the color difference is calculated as follows:

[0146] ;

[0147] For pixel 2, the color difference is calculated as follows:

[0148] ;

[0149] For pixel 3, the color difference is calculated as follows:

[0150] ;

[0151] Calculate the color volatility index of region y:

[0152] ;

[0153] The color volatility index can intuitively reflect the color stability of the region. The higher the value, the greater the color volatility. Finally, the volatility value is used as the basis for simulating the terrain surface features for further 3D reconstruction. During the reconstruction process, regions with greater volatility represent greater terrain changes, while regions with less volatility represent relatively flat terrain.

[0154] The above is only a preferred embodiment of the present invention, and it does not limit the present invention in other forms. Any person skilled in the art may use the disclosed technical content to make changes or modifications into equivalent embodiments with equivalent changes and apply them to other fields. However, as long as it does not depart from the technical solution content of the present invention, any simple modification, equivalent change, and modification made to the above embodiments based on the technical essence of the present invention still fall within the protection scope of the technical solution of the present invention.

Claims

1. A method for terrain recognition and three-dimensional reconstruction based on remote sensing mapping images, characterized in that: The following steps are involved: S1: Based on the surface reflectance spectrum data in the remote sensing mapping image, the spectral intensity under the difference band is recorded, and the band with abnormal spectral intensity is identified according to the band number and light intensity, and summarized as the preferred band index; S2: Based on the preferred band index, adjust the spectral filter configuration, read the current configuration parameters of the hyperspectral sensor, test the adjustable range of the spectral response segment, match the main response segment of each band, retain the intersection segment with the preferred band index, and reset the cutoff band and center band of the spectral filter to generate spectral configuration data; S3: extracting key terrain data from the slope information layer, the curvature grid layer and the elevation gradient layer according to the spectral configuration data, calculating the change in the ratio of slope to curvature, identifying the spatial continuity segments, numbering the segments and mapping their geometric features, and obtaining a terrain feature geometry map; S4: According to the abnormal change area in the terrain feature geometry map, analyze the weak response band associated with the spectral configuration data, record the center position of the area, and pair each position with the band number to integrate into a position band comparison table.

2. The method for terrain recognition and three-dimensional reconstruction based on remote sensing mapping images according to claim 1 is characterized in that: The preferred band index includes spectral intensity anomaly index, spectral adjustment response, and band feature identification; the spectral configuration data includes band adjustment range, filter characteristics, and spectral response configuration; the terrain feature geometry map includes terrain continuity identification results, geometric feature mapping, and spatial change analysis results; the position band comparison table includes center position mark, band response matching results, and terrain anomaly index.

3. The method for terrain recognition and three-dimensional reconstruction based on remote sensing mapping images according to claim 1, characterized in that: The steps for obtaining the preferred band index are specifically as follows: S111: Based on the surface reflectance spectrum data in the remote sensing mapping image, according to the band number and the spectrum intensity, the light intensity and the spectrum intensity difference of each band are identified, and a mapping relationship between the band number and the spectrum intensity is constructed to form a spectrum data mapping structure; S112: using the spectral data mapping structure, calculating the absolute difference between the spectral intensity of each band and the light intensity, identifying the bands exceeding the fluctuation threshold, and determining the band number interval with abnormalities, to obtain the abnormal band interval; S113: Analyze the band numbers in the abnormal band interval, compare the difference between the spectral intensity and the light intensity, classify the band numbers according to the frequency and difference of the band numbers, and summarize them into the preferred band index.

4. The method for terrain recognition and three-dimensional reconstruction based on remote sensing mapping images according to claim 1, characterized in that: The steps for obtaining the spectrum configuration data are specifically as follows: S211: Based on the preferred band index, adjust the spectral filter configuration, align the current center band, cutoff band and response segment boundary parameters of the spectral filter, and select the band segment matching the preferred band by comparing the center and width to obtain the band matching segment; S212: calling the band matching section to refine the response boundary range of each adjustable band, using the formula: ； Get the amount of band overlap , sort the bands, filter the bands with high overlap, and get the band overlap list, where: Indicates Band and The overlap width of the matching segments, It is Band and Alignment factor for matched segments, Indicates Band and The bandwidth interaction characteristics of the matching segments, Indicates Band and Compatibility parameters for matching segments, Indicates the number of matching segments. Indicates the total number of bands; S213: According to the band overlap list, the central band and the cutoff band of each intersection segment are detected, the intersection segment with the preferred band index is retained, and the central band and the cutoff band of the spectral filter are reset to generate spectral configuration data.

5. The method for terrain recognition and three-dimensional reconstruction based on remote sensing mapping images according to claim 1, characterized in that: The steps for obtaining the terrain feature geometric map are specifically as follows: S311: According to the spectral configuration data, the remote sensing mapping image is called, the terrain slope is extracted from the slope information layer, and the terrain slope is matched with the surface curvature extracted from the curvature grid layer, the ratio of the terrain slope to the surface curvature is calculated, and the spatial points that meet the conditions are screened to obtain the slope curvature ratio screening set; S312: analyzing the elevation change of each point according to the slope curvature ratio screening set and combining the elevation gradient data, comparing the consistency of the slope curvature ratio and the elevation change, identifying and numbering the continuous spatial segments, and obtaining the spatial continuity segment number; S313: Based on the spatial continuity segment number, the formula is used: ； Calculate the spatial morphological geometric expression value within the segment , and map its geometric features to obtain the terrain feature geometry map, where Representative In the section The slope change of the point, Representative In the section The elevation fluctuation of the point, Representative In the section The curvature trend of the point, Representative The boundary area of ​​the segment, Representative The edge connectivity of the segment, Representative The mean of the internal distribution density of the segments, Indicates The number of points contained in each segment.

6. The method for terrain recognition and three-dimensional reconstruction based on remote sensing mapping images according to claim 1, characterized in that: The steps for obtaining the position band comparison table are specifically as follows: S411: According to the abnormal change area in the terrain feature geometric map, a boundary undulation section, a slope mutation section and a texture fracture section are selected, image coordinate data and spatial position information of the corresponding sections are collected, spectral configuration data is called, image coordinates and band reflectance sequences are matched, and spatial reflectance corresponding data are obtained; S412: Based on the spatial reflectivity corresponding data, analyze the band reflectivity sequence, calculate the absolute difference between the band reflectivity at each group of coordinate points and the regional average reflectivity, select the band numbers with small reflectivity differences, and record the image coordinate points corresponding to them to obtain the positions of the bands with weak responses; S413: Based on the weak response band position, a pairing relationship between the coordinate point and the band number is constructed and integrated into a position band comparison table.

7. The method for terrain recognition and three-dimensional reconstruction based on remote sensing mapping images according to claim 1, characterized in that: The steps also include: S5: Based on the position band comparison table, extract the pixel set of the associated area, analyze the natural color layer and the surface texture layer of the area, calculate the RGB color mean and local fluctuation, compare the texture color samples, and simulate the terrain surface characteristics to perform three-dimensional reconstruction of the terrain and obtain a three-dimensional terrain model; The three-dimensional terrain model includes three-dimensional coordinate mapping, terrain surface simulation information, and structural hierarchical display.

8. The method for terrain recognition and three-dimensional reconstruction based on remote sensing mapping images according to claim 7 is characterized in that: The steps of obtaining the three-dimensional terrain model are specifically as follows: Based on the position band comparison table, a pixel set of the associated area is extracted, RGB channel data of the natural color layer and the surface texture layer are read, and the average values ​​of the RGB channels are calculated respectively to obtain a regional channel color benchmark; Based on the regional channel color benchmark, the color changes between adjacent pixels in each region are analyzed using the formula: ； Calculating Color Fluctuation , simulate the terrain surface features, reconstruct the terrain 3D, and obtain the terrain 3D model, where: It is The number of pixels in the region, , and It is In the region Red, green, and blue channel data for each pixel, , and It is The average color data of the RGB channels of the region.

Images