Method and system for improving geometric accuracy correction of historical geological maps

Abstract

The present invention discloses a method and system for improving the geometric accuracy correction of historical geological maps. The present invention obtains high-resolution panchromatic images and medium-resolution multispectral images corresponding to historical geological maps, performs geometric correction on the high-resolution panchromatic images with standard remote sensing images as a reference, and then performs geometric correction on the medium-resolution multispectral images with the corrected high-resolution panchromatic images as the reference images; finally, histogram matching and fusion are performed on the single-band intensity images of the high-resolution panchromatic images and the medium-resolution multispectral images to generate historical surface reference color images or further three-dimensional analysis for geometric correction and patch content correction of historical geological maps. The present invention effectively solves the difficult problem of high-precision geometric correction and patch content verification and correction of historical geological maps due to their long age and large surface changes by reconstructing high-precision historical reference scenes, thereby obtaining high-quality historical geological maps, and can perform spatiotemporal change analysis and application with current geological maps.

Description

Technical Field

[0001] The present invention belongs to the field of cartography and remote sensing technology, and in particular relates to a method and system for improving the geometric accuracy correction of historical geological maps. Background Art

[0002] Historical geography maps refer to maps from quite a long time ago (e.g., from the 1970s and 1980s). These maps include soil maps, land use maps, vegetation maps, geological maps, and geomorphological maps. Limited by the data acquisition and mapping methods and techniques of the time, these maps often contain a series of inevitable errors in their spatial location and content. Across agriculture, forestry, environment, ecology, land resources, and water conservancy, there is a widespread demand for the correction, revision, enhancement, and improvement of historical geography maps.

[0003] Take soil maps, for example. Numerous countries around the world have conducted soil surveys, and their most important output is soil maps. Soil maps reflect the soil types and spatial distribution of each region and are the most direct indicator of the quantity and quality of soil resources. For example, during China's Second National Soil Survey, conducted between 1979 and 1986, most townships completed 1:10,000 soil maps, while counties (districts) completed 1:50,000 soil maps. In 2021, the Third National Soil Survey was conducted. The soil maps required for the Third National Soil Survey were refined and updated based on the soil maps obtained from the Second National Soil Survey (referred to as the Second National Soil Survey soil maps). The Third National Soil Survey also analyzed and captured changes in the quantity and quality of soil resources over the past 40 years. Therefore, high-quality Second National Soil Survey soil maps are essential for comparing these soil changes. The general process for producing Second National Soil Survey soil maps for counties, cities, and districts nationwide involves field surveys using topographic maps (or aerial photographs where available), soil profiles, and laboratory analysis and testing data to determine soil types and delineate soil type boundaries. If aerial photographs are used as a base map, stereoscopic observation of adjacent pairs of aerial photographs can be performed under a stereoscope, which helps to better judge and determine soil boundaries. Using aerial photographs as a base map provides a realistic image of the study area, greatly accelerating the speed and quality of soil surveys and offering a significant improvement over methods using topographic maps as a base map. However, using aerial photographs as a base map also brings some new problems, mainly the following:

[0004] 1) Aerial photograph distortion: Aerial photographs are centrally projected, subject to projection errors, tilt errors, and significant edge distortion. Different parts have varying scales. In areas with large terrain variations, elevation differences lead to scale variations, further exacerbating distortion. In areas with large altitude differences, the scales of different elevation areas are inconsistent, leading to distortion.

[0005] 2) Transfer error: A single aerial photograph can only cover a limited area. For a typical county in my country, tens of thousands of aerial photographs are required to fully cover it. Distortion varies from one photograph to another, and even between adjacent photographs.

[0006] 3) Map edge splicing issues: Distortion between different aerial photographs and different map sheets, as well as expansion and contraction of map media (such as paper), lead to misalignment and discontinuity of map patches at the splicing points.

[0007] Therefore, due to the limitations of the aerial photograph soil survey and mapping technology at the time and its related overall soil survey and mapping technology, the final soil map had obvious multiple spatial position errors in the accuracy of the soil boundary positions.

[0008] Over time, especially in recent years, with the rapid economic development, the drastic land use / cover changes have made it extremely difficult to directly use modern high-precision images or topographic maps to align and correct the second-generation soil map. This is mainly reflected in the following aspects:

[0009] 1) The selection of synonyms is extremely difficult or even impossible: This is a fundamental bottleneck in current geometric correction methods when applied to the Second National Soil Map. The Second National Soil Map is a thematic map that primarily displays soil types, with a limited amount of auxiliary information such as rivers, roads, villages, and town boundaries. Nearly 40 years have passed since the Second National Soil Map was completed. During this period, urban expansion has occurred, road and water systems have changed significantly, farmland has been significantly regularized, and the land surface has undergone significant changes. Directly finding a sufficient number of stable synonyms (i.e., feature points that exist in both periods and have not changed in morphological position) with high accuracy, spatial correspondence, and a relatively uniform distribution between the original Second National Soil Map and modern high-precision imagery / topographic maps is extremely difficult and, in most areas, impossible. Without reliable synonyms, traditional point-based geometric registration and correction methods cannot be effectively implemented.

[0010] 2) Existing methods for correcting soil maps from the Second National Survey are severely inaccurate. Even when using existing technical specifications for correction, improving accuracy remains challenging due to a lack of high-quality synonymous points. Random inspections of soil maps from the Second National Survey in selected counties and cities revealed an average geometric error of 50 meters in soil boundaries, with errors reaching as high as 200 meters in some areas. This level of accuracy falls far short of meeting numerous current and future needs.

[0011] 3) Limitations of other data sources: If only contemporary Landsat / SPOT imagery is used, while these images are synchronized with the Second National Geopark (approximately 1979-1986), their spatial resolution (20-80 meters) is insufficient to support meter-level precision positioning and correction, making it difficult to accurately identify subtle features for correction. If only modern high-resolution imagery is used, the aforementioned temporal and spatial mismatches of features persist, making it impossible to find reliable synonyms for accurate registration of historical maps. If only original aerial photographs are used, obtaining them is extremely difficult, and even if original photographs are available, the processing process is cumbersome, requiring a series of steps including scanning, distortion correction, mosaicking, and remapping. This is costly, complex, and often introduces errors when remapping to topographic maps. In particular, individual aerial photographs cover a limited area, making processing large areas extremely time-consuming and costly.

[0012] Of course, in addition to the Second National Soil Map, historical geological maps such as land use maps, vegetation maps, geological maps, and geomorphological maps of the same period (1970s and 1980s) also have similar technical difficulties in performing high-precision geometric correction and improving map content due to their age, huge surface changes, and lack of synonyms. Summary of the Invention

[0013] The present invention aims to address the technical difficulties in performing high-precision geometric correction and image content enhancement on historical geological maps in the prior art, and to provide a method and system for geometrically correcting and enhancing historical geological maps. By integrating multi-source remote sensing data from historical periods, the present invention constructs a high-precision historical surface benchmark color image as an intermediate transitional image information layer. Based on this, the present invention establishes a correlation between current high-precision standard remote sensing images and early historical geological maps, thereby achieving geometrically accurate correction and image content enhancement on the second-generation soil map.

[0014] The specific technical solutions adopted in the present invention are as follows:

[0015] In a first aspect, the present invention provides a method for improving the geometric accuracy of historical geological maps, comprising:

[0016] S1. According to the spatial and temporal range of the historical geological map to be corrected, obtain the corresponding high-resolution panchromatic image and medium-resolution multispectral image, and preprocess them respectively;

[0017] S2. Using the standard geospatial information as the first reference image, selecting control points on the preprocessed high-resolution panchromatic image and the first reference image, and performing geometric correction on the high-resolution panchromatic image;

[0018] S3, using the geometrically corrected high-resolution panchromatic image as the second reference image, selecting control points on the preprocessed medium-resolution multispectral image and the second reference image, and performing geometric correction on the medium-resolution multispectral image;

[0019] S4. Extracting a multi-band image within the spectral range of the high-resolution panchromatic image from the geometrically corrected medium-resolution multispectral image and converting it into a single-band intensity image, then performing histogram matching on the single-band intensity image and the geometrically corrected high-resolution panchromatic image, and then generating a historical surface reference color image through image fusion.

[0020] S5. Using the historical surface benchmark color image as a two-dimensional reference base map, perform geometric correction and image content correction on the historical geological map.

[0021] As a preferred embodiment of the above-mentioned first aspect, after obtaining the historical surface benchmark color image, the historical surface benchmark color image is used as the surface color source, and the digital elevation model data within the said time and space range is used as the elevation source. A three-dimensional reference base map is generated through three-dimensional visualization technology, and the historical geological map is geometrically corrected and the image content is corrected.

[0022] As a preferred embodiment of the first aspect, the historical geological map is a soil map, a land use map, a vegetation map, a geological map or a landform map.

[0023] As a preferred embodiment of the first aspect, the historical geological map is a two-dimensional soil map, the high-resolution panchromatic image is a Conora satellite image, and the medium-resolution multispectral image is a Landsat MSS image, a Landsat TM image or a SPOT HRV image.

[0024] As a preferred embodiment of the first aspect, the geographic reference data adopts a topographic map with a scale of not less than 1:10,000 or an orthophoto with a spatial resolution of not less than 1 meter.

[0025] As a preferred embodiment of the above-mentioned first aspect, the single-band intensity image is obtained by weighted averaging of the multi-band image within the spectral range, or the first principal component of the multi-band image within the spectral range is adopted, or the first component of the Gram-Schmidt transform of the multi-band image within the spectral range is adopted.

[0026] As a preferred embodiment of the first aspect mentioned above, when performing histogram matching, it is necessary to upsample the single-band intensity image to the same resolution as the high-resolution panchromatic image after geometric correction, and then perform histogram matching on the two so that the cumulative histogram of the high-resolution panchromatic image approaches the upsampled single-band intensity image. Finally, the two types of images after histogram matching are fused using the Gram-Schmidt algorithm to generate a historical surface benchmark color image.

[0027] As a preferred embodiment of the first aspect mentioned above, when geometrically correcting and correcting the map content of the historical geological map based on the two-dimensional reference base map and the three-dimensional reference base map, it is necessary to select the same-name landform points of the same period as the control points on the reference base map and the historical geological map, perform geometric correction on the historical geological map, and then compare with the reference base map to correct and improve the boundaries and attribute contents of the soil maps on the geometrically corrected historical geological map, and finally output the final map after topological relationship checking and correction.

[0028] In a second aspect, the present invention provides a computer-readable storage medium having a computer program stored thereon. When the computer program is executed by a processor, the method for improving the geometric accuracy correction of historical geographical maps as described in any one of the schemes in the first aspect above is implemented.

[0029] In a third aspect, the present invention provides a computer electronic device comprising a memory and a processor;

[0030] The memory is used to store computer programs;

[0031] The processor is used to implement the method for improving the geometric accuracy correction of historical geographical maps as described in any one of the solutions of the first aspect when executing the computer program.

[0032] Compared with the prior art, the present invention has the following beneficial effects:

[0033] 1) Solving the problem of historical map correction: This invention reconstructs a high-precision historical benchmark scene by fusing high- and low-resolution remote sensing images from a specific historical period. This effectively solves the long-standing key technical bottleneck of high-precision geometric correction of historical geological maps such as the Second National Soil Map due to the lack of stable homonymous points across time periods due to the long history and large surface changes.

[0034] 2) Significantly Improved Geometric Accuracy: Compared to existing methods, which typically have correction accuracies of only tens of meters or even worse (e.g., 10-200 meters), the proposed method can improve the geometric correction accuracy of historical geological maps, such as the Second National Soil Map, to 2-5 meters. This order-of-magnitude improvement in accuracy lays a solid foundation for subsequent refined applications.

[0035] 3) Significantly improve the accuracy and detail of image content interpretation: This invention not only solves positioning problems but also significantly enhances the content quality of historical geospatial maps, such as soil maps. Using high-precision, information-rich historical color imagery (and 3D views), existing image boundaries can be precisely adjusted, errors corrected, and omissions addressed.

[0036] 4) High-efficiency processing: Compared with processing raw aerial photographs, this method utilizes satellite imagery with large single-view coverage (e.g., the Corona image covers approximately 700 square kilometers). Combining modern geographic information system methods and cloud computing technologies (such as the GEE platform), this method achieves high processing efficiency and is very suitable for large-scale applications.

[0037] 5) Providing Unique Historical Benchmark Data: The high-resolution (meter- and sub-meter-level), multispectral historical surface benchmark color imagery generated by this invention is an extremely valuable data product. It restores the true color landscape of the surface during the historical period corresponding to historical geological maps, filling the gap in high-precision color surface data for that period. The unique combination of Corona and Landsat / SPOT data is key to achieving this goal: Corona provides the necessary meter- and sub-meter-level geometric accuracy framework and micro-relief details, while Landsat / SPOT provides critical color and multispectral information to enhance the interpretability of features, particularly soil-related features. While the human eye can roughly distinguish 10 grayscale levels, this invention can distinguish features with tens of thousands of different color levels. Therefore, this combination is not a simple overlay, but rather an organic fusion, forming a high-quality, information-rich historical feature benchmark.

[0038] 6) Three-dimensional display: The present invention combines color images with meter-level and sub-meter-level geometric accuracy framework and micro-topographic details with DEM (digital terrain elevation model data) to perform three-dimensional display. This is very helpful for interpreting and judging specific soil types and plays an important role in the verification, correction and adjustment of the patch types of the second national soil map.

[0039] 7) Wide range of applications: This method is not only applicable to the second national soil map, but also to the correction, restoration and improvement of other historical geological maps of the same period, such as land use maps, vegetation maps, geological maps, geomorphological maps, etc., or directly used to produce high-precision thematic maps of the period. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 A schematic diagram of the steps for improving the geometric accuracy of historical geography maps;

[0041] Figure 2 It is a schematic diagram of the structure of computer electronic equipment;

[0042] Figure 3 The two-dimensional reference base map and the three-dimensional reference base map generated based on the original Conora satellite image and Landsat TM multispectral image in the embodiment of the present invention;

[0043] Figure 4This is the effect of superimposing a two-dimensional soil map on a two-dimensional reference base map in mountainous and plain areas in an embodiment of the present invention;

[0044] Figure 5 This is a comparison of the soil maps of the exemplary area before and after correction in an embodiment of the present invention. DETAILED DESCRIPTION

[0045] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the specific embodiments of the present invention are described in detail below with reference to the accompanying drawings. In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways than those described herein, and those skilled in the art can make similar improvements without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below. The technical features in the various embodiments of the present invention can be combined accordingly without conflicting with each other.

[0046] In the description of the present invention, it should be understood that the terms "first" and "second" are used solely for descriptive purposes and are not to be construed as indicating or implying relative importance or implicitly specifying the number of technical features being described. Therefore, features defined as "first" or "second" may explicitly or implicitly include at least one of such features.

[0047] The present invention provides a method for geometrically accurately correcting and improving historical geological maps. The method is based on creating a historical surface benchmark color image that fuses multi-source remote sensing data of a specific historical period as an intermediate transition image information layer, thereby accurately correcting historical geological maps and improving the high-quality correction of map content.

[0048] It should be noted that the historical geological maps applicable to the method of the present invention may be early soil maps, land use maps, vegetation maps, geological maps or geomorphological maps, especially geological maps collected around the 1980s or earlier, such as the soil maps obtained by China's second soil survey (i.e., the second soil survey soil map). These geological maps are limited by the technical conditions at the time and often have large errors.

[0049] In a preferred embodiment of the present invention, Figure 1 As shown, the geometric accuracy correction and improvement method of the above historical geological map specifically includes:

[0050] S1. According to the spatial and temporal range corresponding to the historical geological map to be corrected, obtain the corresponding high-resolution panchromatic image and medium-resolution multispectral image, and preprocess them respectively.

[0051] It should be noted that the terms "high resolution" and "medium resolution" in this application refer to relative spatial resolutions and do not limit absolute spatial resolution values. Theoretically, the higher the spatial resolution of high-resolution panchromatic images, the better, and the higher the spatial resolution of medium-resolution multispectral images, the better. However, due to the limitations of remote sensing images available during the historical period of historical geological maps, panchromatic and multispectral images of arbitrary spatial resolutions cannot be obtained. Therefore, it is recommended to select remote sensing images with the highest resolution possible within the available range.

[0052] It should be noted that the preprocessing of high-resolution panchromatic and medium-resolution multispectral images here can be implemented based on the data characteristics of remote sensing images and the conventional preprocessing steps for remote sensing images. Common preprocessing operations include contrast stretching, histogram equalization, radiometric calibration, atmospheric correction, cloud removal, striping noise removal, and bad line processing. Different preprocessing operations can be combined according to the needs of different remote sensing images.

[0053] It should be noted that the temporal and spatial ranges corresponding to the high-resolution panchromatic and medium-resolution multispectral images in the present invention should be as close as possible to the temporal and spatial ranges corresponding to the historical geological map to be corrected. In theory, it is optimal if the temporal and spatial ranges of the three are completely consistent. However, due to the limitations of remote sensing data resources, it is not completely guaranteed that the corresponding remote sensing images can be obtained at any time, and image data may be missing at some locations in space. However, since the distribution of surface features (soil distribution, land use type distribution, landform distribution, etc.) does not change in real time, even if there are certain deviations in the temporal range, it will not have a significant impact. Therefore, the temporal ranges of the high-resolution panchromatic and medium-resolution multispectral images should be as close as possible to the temporal range corresponding to the historical geological map to be corrected. They do not need to be exactly the same. However, the spatial ranges of the high-resolution panchromatic and medium-resolution multispectral images should be as consistent as possible with the spatial range corresponding to the historical geological map to be corrected. Only when the remote sensing image is incomplete should local images be considered for correction.

[0054] S2. Using the standard geographic spatial information as the first reference image, selecting control points on the preprocessed high-resolution panchromatic image and the first reference image, and performing geometric correction on the high-resolution panchromatic image.

[0055] It should be noted that the above-mentioned geographic reference data is used as a reference image for geometric correction of the high-resolution panchromatic image, so its own accuracy needs to meet the correction requirements. It should also be noted that the corresponding time and space ranges of the above-mentioned high-resolution panchromatic image and medium-resolution multispectral image should be as close as possible to the time and space ranges corresponding to the historical geological map to be corrected, but the above-mentioned geographic reference data is not subject to this restriction. The latest accurate data can be used. In theory, the higher the accuracy, the better. In an embodiment of the present invention, the above-mentioned standard geographic spatial information can be a topographic map with a scale of not less than 1:10,000, or a remote sensing orthophoto with a spatial resolution of not less than 1 meter, such as an ESRI global image base map with no offset and a resolution of not less than 1 meter.

[0056] It should be noted that control points (GCPs) refer to paired ground features (also known as keypoints) collected on two images. This step requires selecting keypoints on both the preprocessed medium-resolution multispectral image and the second reference image. Selected control points should be evenly distributed throughout the image, with preference given to features with clear outlines, stable positions, and easy-to-identify locations. These include road intersections, bridges, corners of large, independent buildings, river inflection points, and reservoir dam features. These features existed at the time of the historical map and have remained relatively unchanged since the time of the reference image. Control point density should vary depending on terrain complexity. For example, in gently sloping plains, a control point density of at least 0.5-1 per square kilometer is recommended. In hilly, mountainous areas with significant terrain undulations (e.g., areas with elevation differences exceeding 50 meters and slopes greater than 15 degrees), the density should be increased to 1.5-3 per square kilometer. Therefore, the specific control point selection can be determined based on the actual impact and is not a limitation.

[0057] Furthermore, the geometric correction of the image based on the selected control points in the present invention can be achieved using existing geometric correction algorithms, such as the cubic spline algorithm and the rubbersheet algorithm. However, different images being corrected have different geometric offset characteristics, so the appropriate geometric correction algorithm must be selected based on the image's data. If the high-resolution panchromatic image is a Corona image, the cubic spline algorithm is preferably used for geometric correction and registration to a first reference image (such as the aforementioned ESRI Global Image Base Map).

[0058] S3. Use the geometrically corrected high-resolution panchromatic image as the second reference image, select control points on the preprocessed medium-resolution multispectral image and the second reference image, and perform geometric correction on the medium-resolution multispectral image.

[0059] Similarly, the control points in step S3 above must be paired and collected on the preprocessed medium-resolution multispectral image and the second reference image. The specific selection principles can also refer to the principles described above. The geometric correction algorithm used in this step can also be selected based on the data conditions of the medium-resolution multispectral image itself. If the medium-resolution multispectral image is a Landsat or SPOT image, it is also recommended to use the cubic spline interpolation function (Spline) algorithm for geometric correction and align it to the geometrically corrected high-resolution panchromatic image.

[0060] S4. According to the spectral range of the high-resolution panchromatic image, a multi-band image within the spectral range is extracted from the geometrically corrected medium-resolution multispectral image and converted into a single-band intensity image. The single-band intensity image is then histogram-matched with the geometrically corrected high-resolution panchromatic image, and then a historical surface benchmark color image is generated through image fusion.

[0061] In an embodiment of the present invention, the single-band intensity image is obtained by weighted averaging of multi-band images within the spectral range of the high-resolution panchromatic image (the weights may also be set to be the same, i.e., performing arithmetic averaging), or by using the first principal component (PC1) of the multi-band image, or by using the first component (GS1) of the Gram-Schmidt (GS) transformation of the multi-band image.

[0062] In an embodiment of the present invention, when performing histogram matching on a single-band intensity image and a geometrically corrected high-resolution panchromatic image, the specific method can be implemented as follows: first, the single-band intensity image is upsampled (using methods such as bilinear interpolation and cubic convolution interpolation) to the same resolution as the geometrically corrected high-resolution panchromatic image, and then the two are histogram matched so that the cumulative histogram of the high-resolution panchromatic image approaches the upsampled single-band intensity image. Finally, the single-band intensity image and the high-resolution panchromatic image after histogram matching are fused using the Gram-Schmidt (GS) algorithm to generate a historical surface benchmark color image.

[0063] There are many types of image fusion in the present invention. The above-mentioned Gram-Schmidt (GS) algorithm is particularly suitable for historical soil maps, such as the second-generation soil map.

[0064] S5. Use the historical surface benchmark color image as a two-dimensional reference base map to perform geometric correction and image content correction on the above historical geological maps.

[0065] It should be noted that geometric correction and patch content correction of the aforementioned historical geography maps are achieved by interpreting information displayed on a reference basemap. Because the historical surface reference color imagery utilizes contemporaneous high-resolution panchromatic and medium-resolution multispectral imagery to restore the true color landscape of the surface during that period, it can fill the gap in high-precision color surface data for that period. Therefore, through visual interpretation or automated interpretation using algorithms or neural network technologies, it is possible to accurately extract landforms and features from the historical surface reference color imagery. This historical surface reference color imagery not only ensures meter-level geometric positioning accuracy but also provides rich feature information (color and texture). In practice, historical geography maps can be overlaid on the historical surface reference color imagery basemap, allowing for intuitive verification and correction of information within the historical geography maps by comparing them against the historical surface reference color imagery. This greatly facilitates and improves the efficiency and quality of precisely adjusting patch boundaries, supplementing missed patches, and accurately verifying and correcting patch attributes.

[0066] In addition, the historical surface reference color image obtained in step S4 is two-dimensional. It can be used as a two-dimensional reference base map to perform geometric correction and patch content correction on historical geological maps, but it is not very intuitive. Therefore, in other embodiments of the present invention, after obtaining the historical surface reference color image, the historical surface reference color image can be used as the surface color source, and the digital elevation model (DEM) data within the temporal and spatial range corresponding to the historical geological map can be used as the elevation source to generate a three-dimensional reference base map through three-dimensional visualization technology to perform geometric correction and patch content correction on the historical geological map. In theory, the higher the accuracy of DEM data, the better. Insufficient DEM accuracy may affect the recognition of micro-topographic features and the accuracy of terrain-based boundary delineation. In this embodiment of the present invention, it is recommended to use DEM data with a spatial resolution better than 30 meters and high elevation accuracy, such as the ALOS PALSAR 12.5m DEM, which is available from public channels. Of course, if DEM data with higher resolution is available, the effect will be even better.

[0067] Therefore, the present invention combines historical surface benchmark color images with meter-level and sub-meter-level geometric accuracy frameworks and micro-topographic details with DEM data, which can display the surface color topography in three dimensions and superimpose historical geological maps on this three-dimensional reference base map. Compared with the two-dimensional reference base map, it is more conducive to the interpretation and judgment of specific land features in most scenarios, and plays an important role in the verification, correction and adjustment of the map type of historical geological maps.

[0068] It should be noted that the above-mentioned two-dimensional reference base map and three-dimensional reference base map in the present invention can be selectively used for geometric correction and patch content correction, or they can be used in combination in different correction links, and the specific usage is not limited. When the historical geological map is geometrically corrected and the patch content is corrected based on the two-dimensional reference base map and the three-dimensional reference base map, it is necessary to select the same-name feature points of the same period as the control points on the reference base map and the historical geological map, perform geometric correction on the historical geological map, and then compare with the reference base map to correct and improve the boundaries and attribute content of the soil patches on the geometrically corrected historical geological map (this can be done manually or through automatic algorithms or neural network models), and finally output the final map after topological relationship inspection and correction. The content of the topological relationship inspection includes at least: no overlap (Overlap) between patches, no gaps (Gap), closed boundaries, and integrity of patch attributes. If errors are found in the topological relationship inspection, they need to be repaired accordingly.

[0069] To better understand the specific implementation process and principles of the present invention, the soil map obtained from the Second National Soil Survey (i.e., the Second National Soil Survey soil map) is used as an example historical geological map to demonstrate the specific methods and technical effects of the aforementioned method for improving the geometric accuracy of historical geological maps. For the Second National Soil Survey soil map, the high-resolution panchromatic imagery used is Conora satellite imagery, and the medium-resolution multispectral imagery used is Landsat MSS, Landsat TM, or SPOTHRV imagery. The method includes multiple steps: data acquisition and preprocessing, geometric correction of high-resolution panchromatic imagery, geometric calibration of multispectral imagery, fusion generation of historical surface benchmark color images, fusion of DEM for 3D visualization, geometric correction of the Second National Soil Survey soil map, and content correction and improvement of the Second National Soil Survey soil map. The specific methods for each step are described below.

[0070] 1) Data acquisition and preprocessing: Obtain high spatial resolution (0.6-2 m) panchromatic imagery from the period of the Second National Soil Survey (approximately 1979-1986). In this embodiment, decrypted Conora satellite imagery is used. Obtain medium spatial resolution (20-80 m) and multispectral (3-6 bands) satellite imagery from the same period. In this embodiment, Landsat MSS / TM imagery (30-80 m spatial resolution, 4-6 spectral bands) or SPOT HRV imagery (20 m spatial resolution, 3 bands) is used. Preprocess the various image data obtained. For Corona images, the film scan images need to be quality checked and contrast enhancement and other processing are performed when necessary. For Landsat and SPOT images, in order to reduce the influence of atmosphere and illumination, radiometric calibration and atmospheric correction are required to convert them into apparent reflectance or surface reflectance (FLAASH, QUAC, 6S model or standard algorithms provided by the corresponding data processing system can be used). At the same time, clouds, cloud shadows, stripe noise or bad lines in the images need to be identified and masked to avoid affecting subsequent processing.

[0071] 2) Geometric correction of high-resolution panchromatic images: Using existing high-precision georeferenced data (in this example, topographic maps with a scale greater than 1:10,000 or high-precision orthophotos with a spatial resolution better than 1 meter) as the first reference image, the acquired high-resolution Corona image is accurately geometrically corrected.

[0072] 2.1) Control Point Selection: Control points should be evenly distributed throughout the image. Priority should be given to features with clear outlines, stable positions, and easy-to-identify features. These features include road intersections, bridges, corners of large stand-alone buildings, river inflection points, and reservoir dam features that existed in the early 1980s and have remained relatively unchanged since the reference data period. Control point density should vary depending on the complexity of the terrain. In gently sloping plains, a control point density of at least 0.5-1 per square kilometer is recommended. In hilly and mountainous areas with significant terrain undulations (e.g., areas with elevation differences exceeding 50 meters and slopes greater than 15 degrees), the density should be increased to 1.5-3 per square kilometer. For a single Corona image (covering approximately 700 square kilometers), it is recommended to select at least 200-300 high-quality control points.

[0073] 2.2) Image Correction: Considering the potential for local nonlinear distortion in Corona images, this implementation uses a cubic spline interpolation algorithm, which can handle local deformation, to geometrically correct the preprocessed Corona images and register them to the first reference image. The correction goal in this implementation is to maintain a root mean square error (RMSE) within 2-5 pixels (approximately 2-5 meters for Corona images). After geometric correction, a quasi-orthorectified Corona image with high geometric accuracy is obtained.

[0074] 3) Multispectral Image Geometric Calibration: Using the quasi-orthorectified Corona image from the previous step as the second reference image, perform geometric calibration (registration) on the contemporaneous medium-resolution multispectral imagery (Landsat and SPOT). Using the same calibration algorithm (Spline) as for the high-resolution panchromatic image, select points with the same name that are clearly visible on both the Corona and medium-resolution multispectral images. This calibrates the medium-resolution multispectral image to the second reference image, ensuring precise spatial alignment of the Landsat and SPOT images with the quasi-orthorectified Corona image.

[0075] 4) Fusion Generation of Historical Surface Reference Color Images: Before fusing the geometrically corrected high-resolution panchromatic image (quasi-orthorectified Corona image) from step 2) with the geometrically corrected medium-resolution multispectral image (Landsat / SPOT) from step 3, radiometric normalization / matching preprocessing is required to align the brightness, contrast, and other radiometric characteristics of the PAN image with the intensity information of the MS image. This allows for better preservation of spectral information and reduced color distortion during subsequent fusion. The radiometric normalization method used in this embodiment is histogram matching, and its specific steps may include:

[0076] (a) Generate a target intensity image: Based on the medium-resolution multispectral image (Landsat / SPOT) geometrically corrected in step 3, generate a single-band target intensity image (i.e., a single-band intensity image). This single-band intensity image can be used to represent the overall brightness of the image. Methods for generating a single-band intensity image from a multispectral image include any of the following methods: a1) to a3):

[0077] a1) Based on the spectral range of the high-resolution panchromatic image, extract the multispectral bands within the spectral range of the geometrically corrected medium-resolution multispectral image (for example, the blue, green, red, and near-infrared bands need to be extracted from the Landsat™ image). Then, perform weighted averaging or direct arithmetic averaging on the multispectral bands to calculate a composite brightness image as a single-band intensity image.

[0078] a2) extracting the first principal component (PC1) of the principal component analysis in the multispectral band obtained in a1) as a single-band intensity image;

[0079] a3) Extract the first component (GS1) of the Gram-Schmidt transform of the multispectral band obtained in a1) as a single-band intensity image.

[0080] It should be noted that different medium-resolution multispectral images, such as Landsat and SPOT, require the extraction of single-band intensity images. The single-band intensity image corresponding to any medium-resolution multispectral image can be used for the subsequent histogram matching and image fusion of the present invention.

[0081] (b) Perform histogram matching:

[0082] The single-band intensity image generated in step (a) has a lower resolution than the high-resolution panchromatic image (quasi-orthorectified Corona image), so it needs to be upsampled to the same spatial resolution as the high-resolution panchromatic image through resampling methods (such as bilinear interpolation and cubic convolution interpolation).

[0083] The high-resolution panchromatic image after geometric correction is used as the source image, and the above-sampled single-band intensity image is divided into the target (reference) image. The histogram matching algorithm is executed to adjust the pixel values of the high-resolution panchromatic image so that its cumulative histogram approaches the cumulative histogram of the target image, and the radiometrically matched panchromatic image (PAN_matched) is obtained.

[0084] After completing the radiation normalization / matching preprocessing, the panchromatic image PAN_matched (high resolution) and the geometrically corrected medium-resolution multispectral image (multispectral bands, maintaining their original resolution) are fused through an image fusion algorithm to generate a high-resolution historical surface benchmark color image that has both high spatial resolution and multispectral information and reflects the surface conditions during the second census period.

[0085] The image fusion algorithm used in this invention preferably uses the Gram-Schmidt (GS) transform (or its modified variants). The GS algorithm is chosen because it effectively incorporates the high spatial resolution details of Corona imagery while effectively preserving the original spectral characteristics of Landsat / SPOT imagery, reducing color distortion during the fusion process. This is crucial for subsequent interpretation of soil-related features (such as vegetation and land use types) that relies on spectral information. In contrast, other fusion methods (such as HSV) can sometimes result in significant color distortion and are therefore not preferred in this invention.

[0086] The fusion parameters of the image fusion algorithm are determined as follows: The specific parameters of the Gram-Schmidt fusion (e.g., weight coefficients) should be determined based on the specific implementation of the software used (e.g., ArcGIS Pro, ENVI, or Google Earth Engine). Optimization can generally be achieved through a combination of qualitative visual interpretation of the fusion results (examining spatial clarity, color naturalness, and the presence of obvious artifacts) and quantitative evaluation metrics (see the next step, quality control).

[0087] After completing the image fusion of the single-band intensity image and the geometrically corrected high-resolution panchromatic image using histogram matching and the Gram-Schmidt algorithm, the resulting historical surface reference color image requires image quality control. Specifically, the quality of the fused historical surface reference color image should be evaluated using a combination of qualitative visual inspection (spatial texture clarity, spectral fidelity, color naturalness, and the presence of fusion artifacts) and quantitative metrics. Quantitative metrics include the relative error of global dimensions (ERGAS), spectral angle mapping (SAM), correlation coefficient (CC), and root mean square error (RMSE). The fused image should be compared with the original multispectral image to assess the degree of spectral distortion and with the original panchromatic image to assess the effect of spatial detail injection. Acceptable thresholds for these metrics should be set to ensure that the fused image quality meets the requirements of subsequent applications.

[0088] 5) DEM Fusion for 3D Visualization: The historical surface reference color image generated in the previous step is fused with digital elevation model (DEM) data to create a 3D visualization of the target area in GIS software (such as ArcGIS Pro). When creating the 3D visualization, the historical surface reference color image should be used as the surface color source, and the digital elevation model data as the elevation source. The spatial extent of the DEM data should correspond to the spatial extent of the historical surface reference color image itself. There is no time limit, and the latest DEM data can be used. The spatial resolution and elevation accuracy of the DEM used for 3D visualization will affect the 3D visualization effect and the accuracy of terrain-based interpretation. This example recommends using DEM data with a spatial resolution better than 30 meters (such as the ALOS PALSAR 12.5m DEM or higher) and high elevation accuracy. Insufficient DEM accuracy may affect the identification of micro-relief features and the accuracy of terrain-based boundary delineation. The 3D visualization scene generated using 3D visualization technology can then serve as a 3D reference basemap for subsequent geometric correction and revision of the second-level soil map.

[0089] 6) Geometric Correction of the Secondary General Soil Map: Using the high-precision historical surface benchmark color image generated in Step 4 (or combined with the 3D visualization scene from Step 5) as a precise reference basemap, perform geometric correction on the original Secondary General Soil Map. Select points on the Secondary General Soil Map that have the same name and are clearly identifiable on the historical surface benchmark color image (note: this search is performed on images from the same historical period, not across 40 years). Geometric correction of the Secondary General Soil Map is then performed. The preferred geometric correction algorithm for the Secondary General Soil Map is the Rubbersheet algorithm, which locally adjusts vector data based on control points to accommodate irregular map deformations. The goal is to achieve meter-level geometric accuracy (e.g., 2-5 meters) for the corrected soil map.

[0090] 7) Correction and Improvement of the Second National Soil Map: The geometrically corrected Second National Soil Map from step 6 is overlaid on the historical surface reference color image from step 4 (or the 3D visualization scene from step 5). Soil experts or trained technicians will correct and improve the boundaries and content of the patches. Automated algorithms can also be designed to correct and improve the boundaries and content of patches.

[0091] The revision and improvement of the second-level soil map content must comply with certain revision standards. In the embodiment of the present invention, the following boundary adjustment basis and content update standard requirements can be set:

[0092] Boundary Adjustment Basis: In a 2D or 3D environment, soil patch boundary adjustments are primarily based on: 1) geomorphic unit boundaries (e.g., ridgelines, valley lines, slope / aspect change lines, and terrace edges) as shown on the fused high-resolution historical color imagery; 2) distinct boundaries of vegetation or land use types as reflected by color and texture on the imagery (this correlation with soil type should be determined in conjunction with soil knowledge); and 3) topographic features (e.g., sunny / shady slopes, top / shoulder / mid / foot of slopes) as shown in conjunction with the DEM. The goal is to ensure that the adjusted soil boundaries closely match the restored 1980s landscape.

[0093] Content update: Verify the accuracy of soil type labeling in the original map by comparing it with historical benchmark images and correct any errors; identify and outline soil patches that were omitted from the original map or could not be expressed due to scale limitations; merge or subdivide the original map patches as needed.

[0094] 8) Topological check and final output: Check and correct the topological relationship of the revised and improved second-level soil map.

[0095] Topological relationship rules: Topological checks should at least include ensuring that patches have no overlaps, gaps, closed boundaries, and complete patch attributes. Topological rule validation tools in GIS software (such as ArcGIS Pro) can be used to perform these checks and manually or automatically fix any errors found.

[0096] Final result quality evaluation: The quality evaluation of the completed soil map correction results of the second national census should include: 1) Geometric accuracy: by selecting independent check points (ICPs), measuring coordinates on the historical benchmark color image, and calculating the RMSE of the corrected map boundary position, which should meet the design requirements (such as 2-5 meters); 2) Topological consistency: confirm that there are no topological errors; 3) Content rationality: soil experts will conduct sampling inspections and evaluations on the geographical / geomorphological conformity of the map boundary and the accuracy of the map attribute annotation by comparing the historical benchmark image and the 3D view.

[0097] Final output: Output digital results of the second national soil census that meet standards, undergo high-precision correction and have complete content.

[0098] This demonstrates that the present invention, by fusing high- and low-resolution remote sensing images from a specific historical period to reconstruct a high-precision historical surface benchmark color image, effectively addresses the long-standing technical bottleneck of high-precision geometric correction of historical maps, such as the Second National Soil Map, due to the lack of stable, cross-period synonymous points caused by the age and significant surface changes. This historical surface benchmark color image establishes a link between current high-precision remote sensing imagery and early Second National Soil Maps, enabling the geometric correction, content revision, and improvement of Second National Soil Maps using the rich information provided by the historical surface benchmark color image.

[0099] Of course, although the above embodiment is described using the second-level soil map as a typical example, it is obvious that the method of the present invention is not only applicable to the second-level soil map, but is also applicable to correcting, repairing and improving other historical geological maps of the same period, such as land use maps, vegetation maps, geological maps, landform maps, etc., or directly used to produce high-precision thematic maps of the period, without limitation.

[0100] It should be noted that the method steps shown in S1 to S5 above can essentially be implemented in the form of a computer program.

[0101] Therefore, based on the same inventive concept, Figure 2 As shown, the present invention also provides a computer electronic device corresponding to the method for geometrically correcting and improving historical geographical maps provided in the above embodiment, which includes a memory and a processor;

[0102] The memory is used to store computer programs;

[0103] The processor is configured to implement the aforementioned method for geometrically correcting and improving historical geographical maps when executing the computer program;

[0104] Furthermore, the logic instructions in the aforementioned memory can be implemented as software functional units and, when sold or used as independent products, can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, or the portion that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions for causing a computer device (which can be a personal computer, server, or network device, etc.) to execute all or part of the steps of the method described in each embodiment of the present invention.

[0105] Therefore, based on the same inventive concept, the present invention provides a computer-readable storage medium corresponding to a method for geometrically correcting and improving historical geographical maps, on which a computer program is stored. When the computer program is executed by a processor, the method for geometrically correcting and improving historical geographical maps as described above can be implemented.

[0106] Therefore, based on the same inventive concept, the present invention provides a computer program product, including a computer program / instruction, which, when executed by a processor, can implement the method for improving the geometric accuracy correction of historical geographical maps as described above.

[0107] Specifically, in the computer-readable storage medium of the above three embodiments, the stored computer program is executed by the processor to perform the above steps S1 to S4.

[0108] It is understood that the storage medium may include random access memory (RAM) or non-volatile memory (NVM), such as at least one disk storage device. Furthermore, the storage medium may be any medium capable of storing program code, such as a USB flash drive, a mobile hard drive, a magnetic disk, or an optical disk.

[0109] It is understandable that the above-mentioned processor can be a general-purpose processor, including a central processing unit (CPU), a network processor (NP), etc.; it can also be a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA) or other programmable logic devices, discrete gate or transistor logic devices, and discrete hardware components.

[0110] It should also be noted that those skilled in the art can clearly understand that for the convenience and brevity of description, the specific working process of the system described above can refer to the corresponding process in the aforementioned method embodiment, and will not be repeated here. In the various embodiments provided in this application, the division of steps or modules in the system and method is only a logical function division. In actual implementation, there may be other division methods, for example, multiple modules or steps can be combined or integrated together, and a module or step can also be split.

[0111] In order to better understand the technical effects that can be achieved by the above-mentioned method of the present invention, the specific implementation method of the geometric precision correction and improvement method of the historical geological map in the above-mentioned embodiment, as well as the correction and improvement effects are demonstrated below through an exemplary local spatial area.

[0112] Example

[0113] In this embodiment, the historical geological map is a two-dimensional soil map. The corresponding high-resolution panchromatic image is a Conora satellite image, and the medium-resolution multispectral image is a Landsat MSS / TM image or a SPOT HRV image. The specific steps of the geometric precision correction and improvement method of the two-dimensional soil map are as follows:

[0114] Step 1: Data acquisition and preprocessing

[0115] The third batch of decrypted Conora satellite images with high spatial resolution (0.6-2 meters) from the Second National Soil Survey (1979-1986) were obtained, as well as remote sensing satellite images with medium spatial resolution (20-80 meters) and multispectral bands (3-6 bands) from the same period. The data were preprocessed. The specific data and preprocessing methods are as follows:

[0116] Obtain high-resolution panchromatic imagery: Obtain declassified Corona KH series imagery (such as KH-4B and KH-9, with a spatial resolution of approximately 0.6-2 meters) from the U.S. Geological Survey (USGS) EarthExplorer website. Specifically, select Declassified Data -> Declass 3 in DataSets and download it.

[0117] Obtaining medium-resolution multispectral imagery: Obtain contemporaneous Landsat MSS (approximately 80-meter resolution, 4 bands) or Landsat TM (30-meter resolution, 6 visible / infrared bands) imagery (Landsat Collection 2 Level-1) from USGS EarthExplorer; or obtain SPOT-1 HRV imagery (20-meter multispectral resolution, 3 bands) from the SPOT website. If multiple images are available for the same period, select the one with less cloud cover and better quality.

[0118] Corona image data preprocessing: Visually inspect the downloaded scanned image to confirm that it is free of major scratches or stains. Use image processing software (such as Photoshop, GIMP, or tools built into your GIS software) to perform contrast stretching or histogram equalization to enhance detail.

[0119] Landsat / SPOT image data preprocessing: Landsat / SPOT image data can be processed using professional remote sensing image processing software (such as ENVI, ERDAS IMAGINE, ArcGIS Pro, QGIS) or online platforms (such as Google Earth Engine - GEE). Specific preprocessing includes radiometric calibration, atmospheric correction, and cloud and noise processing, as follows:

[0120] Radiometric calibration: Convert DN values to radiance.

[0121] Atmospheric correction: Use the surface reflectance product algorithm provided by FLAASH, QUAC, 6S model or GEE to convert radiance or apparent reflectance into surface reflectance to eliminate the effects of atmospheric scattering and absorption.

[0122] Cloud and noise processing: Use the image's built-in quality assessment band (QA band) or specialized cloud detection algorithms (such as Fmask) to identify and generate masks for clouds and cloud shadows. Additionally, any existing banding noise or bad lines should be repaired (e.g., using a debanding tool).

[0123] Step 2: High-resolution panchromatic image high-precision geometric correction

[0124] In the embodiment, Conora images of the required location and time are selected and downloaded from USGS, and loaded into professional geospatial commercial software, such as ArcGIS Pro, MapGIS, SuperMap, etc., using an ESRI global image base map with a resolution higher than 1 meter without offset.

[0125] Select control points (GCPs): Manually select points with the same name in the Corona image and reference data, following the control point selection principles detailed above (even distribution, stable and clear features) and density requirements (differentiated by terrain, e.g., 0.5-1 / km² in plains, 1.5-3 / km² in mountainous areas, totaling 200-300+). Record the image coordinates and corresponding geographic coordinates (or reference image coordinates) for each control point.

[0126] Select the calibration model and execute: Select the cubic spline algorithm as the geometric transformation model. Enter the selected GCPs and execute the calibration. Accuracy Assessment: Check the RMSE value after calibration to ensure it meets the target accuracy (e.g., 1-3 meters). After completing the geometric calibration of the calibration model, a quasi-orthorectified Corona image is obtained.

[0127] Step 3: Multispectral Image (Landsat / SPOT) Geometric Calibration

[0128] In the GIS software, geometric calibration is performed using the quasi-orthorectified Corona image corrected in step 2 as the reference image and the multispectral image (Landsat / SPOT) as the image to be calibrated. In this example, Landsat-5 TM images are used as an example image to be calibrated (the calibration process for SPOT is similar). The specific calibration process is as follows:

[0129] Select key points: Select key points between the multispectral image to be calibrated and the quasi-orthorectified Corona image.

[0130] Perform calibration: Using the cubic spline interpolation function (Spline) algorithm as the geometric transformation model, the multispectral image is geometrically calibrated (registered) based on the selected congruent points. This ensures that the calibrated multispectral image and the Corona image are accurately spatially superimposed.

[0131] Step 4: Fusion generation of high-resolution historical surface benchmark color images

[0132] This step combines the spatial details of the high-resolution panchromatic image (PAN, i.e., the Corona image with a spatial resolution of 1 meter) calibrated in Step 2 with the spectral information of the medium-resolution multispectral image (MS, i.e., the Landsat-5 TM image with a spatial resolution of 30 meters) calibrated in Step 3. This creates a historical surface reference color image that combines high spatial resolution (1 meter) with multispectral information (inheriting the bands of Landsat TM). This step can be performed locally using professional remote sensing / GIS software (such as ENVI and ArcGIS Pro) or cloud platforms (such as Google Earth Engine (GEE)). Given the potentially large data volumes of Corona and Landsat imagery, utilizing the computing power of cloud platforms like GEE offers efficiency advantages for processing large areas.

[0133] This embodiment uses the Gram-Schmidt (GS) transform fusion algorithm and performs radiometric normalization (using histogram matching as an example) as a preprocessing step before fusion to optimize the fusion effect, especially color fidelity. The specific fusion process is as follows:

[0134] 1) Radiometric Normalization - Histogram Matching

[0135] (a) Generating a Target Intensity Image: A single-band intensity image representing the overall brightness of the original 30-meter-resolution Landsat-5™ surface reflectance image is calculated. This example uses arithmetic averaging to average Landsat-5™ Bands 1 (blue), 2 (green), 3 (red), and 4 (near-infrared), which overlap well with the Corona spectral range (visible to partial near-infrared). The averaged single-band intensity image, Intensity_MS, is calculated using the formula: Intensity_MS = (Band1_reflectance + Band2_reflectance + Band3_reflectance + Band4_reflectance) / 4, where Band1_reflectance, Band2_reflectance, Band3_reflectance, and Band4_reflectance represent Band 1 (blue), Band 2 (green), Band 3 (red), and Band 4 (near-infrared) in the Landsat-5™ image, respectively. The calculations were performed at the native 30-meter resolution of the Landsat imagery.

[0136] (b) Upsampling the target intensity image: Using the bilinear interpolation method, the 30-meter resolution Intensity_MS image generated in the previous step is resampled to the same spatial resolution as the Corona image (1 meter), thereby obtaining the high-resolution target intensity image Intensity_MS_resampled.

[0137] (c) Perform histogram matching: Using the original high-resolution Corona image (PAN) corrected in step 2 as the source image (Source) and Intensity_MS_resampled as the target / reference image (Target / Reference), perform histogram matching using image processing software (such as the Histogram_Match tool in ENVI or the skimage.exposure.match_histograms function in Python). This process adjusts the pixel values (grayscale) of the original Corona image so that its cumulative frequency histogram closely matches that of Intensity_MS_resampled. The resulting radiometrically matched panchromatic image, Corona_matched, retains the original 1-meter spatial resolution and ground detail of the Corona image, but its overall brightness and contrast distribution has been adjusted to more closely resemble the intensity characteristics of Landsat imagery.

[0138] 2) Gram-Schmidt (GS) fusion:

[0139] Image fusion can be implemented on local GIS / remote sensing software (such as ArcGIS Pro's "Create Pan-sharpenedRaster Dataset" tool selecting the GS method, or the corresponding module of ENVI) or cloud platforms (such as GEE).

[0140] Local software example (ArcGIS Pro): Use the Create Pan-sharpened Raster Dataset tool, input the rectified Corona image (as the panchromatic band) and the calibrated Landsat™ multispectral image (as the multispectral band), and select the Gram-Schmidt fusion method. Adjust the weight parameters as needed (the software typically has default values or provides optimization options).

[0141] Cloud Platform Example (GEE): Upload the corrected / calibrated image to GEE Assets. Write a script to call the pan-sharpening function provided by GEE (you may need to find a function or community script that supports the GS algorithm) and perform the fusion. GEE's strengths lie in its ability to handle large amounts of data and parallel computing.

[0142] Parameter Optimization and Quality Assessment: Visually inspect the fused image for spatial clarity, color naturalness, and the presence of obvious artifacts. Quantitative metrics (such as ERGAS, SAM, CC, and RMSE) are calculated and compared with the original image. Fusion parameters are adjusted (if the software allows) to achieve optimal results. Ensure that the fused image quality meets pre-set standards. Finally, the fused image is output as a high-resolution, multispectral, historical surface benchmark color image.

[0143] Step 5: Fuse DEM for 3D visualization

[0144] In this embodiment, the historical surface reference color image generated in step 4 is loaded into ArcGIS Pro, a GIS software that supports 3D scenes. DEM data corresponding to the area is also loaded (ALOS 12.5m DEM or higher-precision data is recommended, and preprocessing is required to align the image coordinate system). The historical surface reference color image is set as the surface color source for the 3D scene, and the DEM is set as the elevation source for the 3D scene. A local scene (Local Scene) or a global scene (GlobalScene) is created for 3D display. This generates a 3D reference base map using 3D visualization technology, allowing users to freely zoom, rotate, and tilt the viewing angle using the mouse for stereoscopic observation.

[0145] like Figure 3 The figure shows a 2D reference basemap (historical surface benchmark color image) and a 3D reference basemap (3D visualization scene) generated based on original Conora satellite imagery and Landsat™ multispectral imagery. As can be seen, the original Conora (black and white) satellite imagery has high spatial resolution and clear feature details. However, black and white images are central projections without geographic coordinates, resulting in inconsistent scales for different parts and overall image distortion. The Landsat™ multispectral (color) imagery, on the other hand, has a spatial resolution of 30 meters, but features are difficult to discern and exhibit geometric distortion. After fusing the two images according to the method of this embodiment, the resulting historical surface benchmark color image combines high spatial resolution with multispectral information (partial color image), clearly displaying information on different land types and their boundary details. Furthermore, by creating a 3D visualization scene with the aid of DEM information, soil and feature types can be more easily displayed, observed, analyzed, and judged from different angles and directions, significantly facilitating soil map correction and quality improvement.

[0146] Step 6: Geometric correction of the second soil map

[0147] In this embodiment, after overlaying the two-level soil map in the GIS software, the historical reference color image completed in step 4 is used as a precise reference basemap. Using georeferencing tools (such as the Georeferencing toolbar or the SpatialAdjustment tool in ArcGIS Pro), identical features (such as road intersections, river confluences, reservoir corners, and unique landforms that are clearly identifiable in both) are selected between the two-level soil map and the historical reference image. The Rubber Sheet transformation model is selected to perform geometric correction. Accuracy Verification: Independent Check Points (ICPs) are used to assess the correction accuracy, with the goal of keeping the positional error of the two-level soil map's patch boundaries within 2-5 meters.

[0148] like Figure 4 As shown in Figures (a) and b), respectively, illustrate the effects of overlaying a two-dimensional soil map with a two-dimensional historical color image in mountainous and plain areas. As can be seen, the deviation (displacement) of soil type (feature) boundaries in mountainous areas can often exceed 60 meters, and in some areas, reach 120 meters. After overlay, using high-spatial-resolution, multispectral-fused imagery as the basemap facilitates the selection of place names and features, geometric correction of the two-dimensional soil map, and easy adjustment of feature boundaries and verification of feature types. Similarly, the displacement (deviation) of soil maps (feature types) in plain areas is typically 30-40 meters. After overlay, using high-spatial-resolution, multispectral-fused satellite imagery as the basemap facilitates the adjustment of feature boundaries and verification of feature types.

[0149] Step 7: Correction and improvement of the second soil map

[0150] In the GIS software, overlay the corrected soil map from step 6 onto the historical surface color image from step 4 (or the 3D visualization scene from step 5) to revise and improve the content of the soil map. This process can include starting editing, adjusting boundaries and adding new content, and verifying and revising content. The steps are as follows:

[0151] Start Editing: Start an editing session for the Soil Map Vector Data.

[0152] Boundary Adjustment and Addition: Using editing tools, we precisely adjust the boundaries of existing soil patches to align with historical surface landscape characteristics, comparing them to the basemap imagery (and 3D scene) and following the aforementioned boundary adjustment criteria and content update standards. For soil units that were omitted from the original image or were not depicted due to scale limitations, we added new boundaries and patches.

[0153] Content Verification and Correction: Check the soil type attribute annotations for each patch. Correct incorrect annotations by comparing image features with soil knowledge. Merge patches that are too fragmented or have incorrect type interpretations, or subdivide larger patches as needed.

[0154] in addition, Figure 5 The paper further demonstrates the results of a patch comparison between the two soil maps before and after correction in an exemplary scenario. This demonstrates that using high-spatial-resolution, multispectral color satellite imagery, supplemented by DEM information, for three-dimensional observation can effectively calibrate patch boundaries, verify and modify patch content, and complete topological verification of the soil map. Figure 5 The comparison results before and after the correction show that the blue line represents the original soil map, while the red line represents the corrected soil map boundary. This clearly demonstrates that the spatial position boundaries of the features in the soil map after spatial geometric position correction generally coincide with the features. In this scene, the total number of patches in the displayed area is 16. Using the fused color image and terrain digital elevation model display for interpretation, eight new feature boundaries and eight new patches were added. Three of these boundaries correspond to three patches that were omitted during the original mapping scale, and five feature boundaries correspond to five patches that were newly added based on the current resolution. Two patch boundaries were deleted, and one patch was merged. This shows that the spatial position correction, patch content verification, and revisions have significantly improved the quality of the second-level soil map for this area.

[0155] Step 8: Topology check and final output

[0156] In this embodiment, topology rules (such as Must Not Overlap, Must Not Have Gaps, etc.) are further created for the corrected soil layer in the second-level soil map in GIS software (ArcGIS Pro), and all topology errors are found and fixed by running topology validation.

[0157] Final quality evaluation: Conduct a final quality check, including three aspects: geometric accuracy, topological consistency, and content rationality:

[0158] Geometric accuracy: ICPs were used to reconfirm that the RMSE was within the 2-5 meter range.

[0159] Topological consistency: Confirm that there are no topological errors.

[0160] Content plausibility: Soil experts conduct spot checks to assess boundary plausibility and attribute accuracy.

[0161] After passing the final quality inspection, the final results can be output and the final high-precision, complete digital results of the second-level soil map can be saved (such as shapefile, geodatabase feature class and other standard formats).

[0162] In addition, two comparative experiments are designed in this embodiment to verify the superiority of using the high-precision historical reference color image of the present invention to correct and enhance the second-level soil map.

[0163] 1. Experimental area and data

[0164] The longitude and latitude range of the experimental area is: 120.5315-120.5684 east longitude, 28.0844-28.1071 north latitude

[0165] The raw data obtained for the experimental area are as follows:

[0166] 1) Original paper or scanned soil map of the experimental area during the Second National Soil Survey; 2) Corona KH series high-resolution panchromatic imagery (spatial resolution approximately 1-2 meters) covering the area in November 1983 during the Second National Soil Survey; 3) Landsat TM multispectral imagery (spatial resolution 30 meters) covering the area during the same period; 4) Modern high-precision reference data ESRI basemap for initial geometric correction; 5) Digital Elevation Model (DEM), ALOS PALSAR 12.5m DEM.

[0167] The evaluation indicators used in the experiment are as follows:

[0168] Geometric correction accuracy: Using independent check points (ICPs), the root mean square error (RMSE) of the corrected two-level soil map patch boundaries relative to the "true" historical surface position (using the high-precision historical color reference image constructed by the method of the present invention or the high-precision Corona image based on it as a reference) is calculated.

[0169] The content improvement effect is evaluated using two dimensions: qualitative and quantitative:

[0170] Qualitative assessment: A soil science or remote sensing interpretation expert, with the assistance of their respective reference images (or 3D views), will evaluate the consistency of the patch boundaries of the calibrated soil map with the landforms / ground features, the internal homogeneity of the patches, and the reliability of the patch attribute interpretation (excellent, good, fair, or poor).

[0171] Quantitative evaluation (reference): During the statistical content revision process, the degree of patch boundary adjustment, the number and area of newly added / deleted / merged patches, the number of attribute revisions, etc.

[0172] 2. Comparative Experimental Design

[0173] The following two comparative experiments were set up. Except for the core reference images, other processing steps (such as initial image preprocessing, the transformation model used for birefringence geometric correction - rubber stretch transformation, the operation process of content correction, topology check, etc.) were kept as consistent as possible with the present embodiment.

[0174] Comparative Experiment 1: Using Only Historical Medium-Resolution Multispectral Imagery as a Baseline

[0175] (1) Acquire and preprocess Landsat TM images of the same period (radiometric calibration, atmospheric correction, etc.).

[0176] (2) Geometric correction of Landsat TM images using modern high-precision reference data.

[0177] (3) The calibrated Landsat TM image (30-meter resolution) is used as the base map.

[0178] (4) Select the same-name points on the base map and perform geometric correction (rubber stretch transformation) on the original two-dimensional soil map.

[0179] (5) Correct and improve the content of the second soil map based on the reference base map (which can be combined with DEM to construct a three-dimensional view).

[0180] (6) Conduct topology inspection, output results, and conduct evaluation.

[0181] Comparative Experiment 2: Using Only Historical High-Resolution Panchromatic Images as a Baseline

[0182] (1) Acquire and preprocess the Corona images of the same period (contrast enhancement, etc.).

[0183] (2) Use modern high-precision reference data to perform high-precision geometric correction on the Corona image (e.g., Spline transformation, the steps are the same as step 2 of the present invention).

[0184] (3) The calibrated Corona image (1-2 m resolution, black and white) is used as the base map.

[0185] (4) Select the same-name points on the base map and perform geometric correction (rubber stretch transformation) on the original two-dimensional soil map.

[0186] (5) Correct and improve the content of the second soil map based on the reference base map (which can be combined with DEM to construct a three-dimensional view).

[0187] (6) Conduct topology inspection, output results, and conduct evaluation.

[0188] 3. The experimental results and analysis are as follows:

[0189] Table 1 Experimental results of different experimental groups

[0190]

[0191] From the experimental results in Table 1, it can be seen that the base map obtained in the comparative experiment 1 has a low resolution and cannot provide fine boundary information of land objects, resulting in insufficient geometric positioning accuracy, and it is difficult to accurately identify and outline small patches and accurately adjust the boundaries; the base map obtained in the comparative experiment 2 has a high spatial resolution and good geometric positioning accuracy. However, the lack of color and spectral information makes it difficult to interpret surface differences caused by vegetation, water bodies, and different soil types, which seriously affects the accurate verification and correction of the content (especially attributes) of the patches. The base map obtained in this embodiment (i.e., the historical surface benchmark color image) combines high spatial resolution and multispectral information, which not only ensures meter-level geometric positioning accuracy, but also provides rich land object information (color, texture), greatly facilitating and improving the efficiency and quality of precise adjustment of patch boundaries, supplementation of missed patches, and accurate verification and correction of patch attributes.

[0192] The results of the above comparative examples clearly demonstrate that using only historical medium-resolution multispectral imagery cannot meet the accuracy requirements for meter-level geometric correction of the Second National Soil Map, and the content correction effect is poor. Using only historical high-resolution panchromatic imagery, while achieving high geometric positioning accuracy, significantly reduces the accuracy and reliability of content correction (particularly attribute interpretation) due to the lack of spectral information. The technical solution proposed in this invention, which constructs and utilizes fused high-precision historical color reference imagery, is the only effective approach that can simultaneously achieve meter-level geometric correction accuracy and high-quality content correction and improvement. By integrating the advantages of two historical data sources, it overcomes the limitations of a single data source and holds the key to solving the challenges of accurately correcting and improving the quality of Second National Soil Maps.

[0193] Statistics from experimental results across a large number of regions show that, by applying the method of the present invention, this embodiment successfully reduced the average geometric position error of the regional soil map from approximately 50 meters to 2-5 meters. In terms of content, by verifying and correcting it against high-precision historical benchmark imagery, the number of image patches doubled from 8 to 16, significantly improving the accuracy and detail of the soil map. This demonstrates the significant effectiveness of the method of the present invention in correcting and improving the quality of historical geological maps, such as the soil map of the Second National Geographic Survey.

[0194] The embodiments described above are merely some preferred implementations of the present invention and are not intended to limit the present invention. Persons skilled in the art may make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, any technical solution obtained by equivalent substitution or equivalent transformation falls within the scope of protection of the present invention.

Claims

1. A method for improving the geometric accuracy of historical geography maps, characterized by: include: S1. According to the spatial and temporal range of the historical geological map to be corrected, obtain the corresponding high-resolution panchromatic image and medium-resolution multispectral image, and preprocess them respectively; S2. Using the standard geospatial information as the first reference image, selecting control points on the preprocessed high-resolution panchromatic image and the first reference image, and performing geometric correction on the high-resolution panchromatic image; S3, using the geometrically corrected high-resolution panchromatic image as the second reference image, selecting control points on the preprocessed medium-resolution multispectral image and the second reference image, and performing geometric correction on the medium-resolution multispectral image; S4. Extracting a multi-band image within the spectral range of the high-resolution panchromatic image from the geometrically corrected medium-resolution multispectral image and converting it into a single-band intensity image, then performing histogram matching on the single-band intensity image and the geometrically corrected high-resolution panchromatic image, and then generating a historical surface reference color image through image fusion. S5. Using the historical surface benchmark color image as a two-dimensional reference base map, perform geometric correction and image content correction on the historical geological map.

2. The method for improving geometric accuracy of historical geography maps according to claim 1, characterized in that: After obtaining the historical surface benchmark color image, the historical surface benchmark color image is used as the surface color source, and the digital elevation model data within the time and space range is used as the elevation source. A three-dimensional reference base map is generated through three-dimensional visualization technology, and the historical geological map is geometrically corrected and the map content is corrected.

3. The method for improving geometric accuracy of historical geography maps according to claim 1, characterized in that: The historical geological map is a soil map, a land use map, a vegetation map, a geological map or a landform map.

4. The method for improving geometric accuracy of historical geography maps according to claim 1, characterized in that: The historical geological map is a two-dimensional soil map, the high-resolution panchromatic image is a Conora satellite image, and the medium-resolution multispectral image is a Landsat MSS image, a Landsat TM image, or a SPOT HRV image.

5. The method for improving geometric accuracy of historical geography maps according to claim 1, characterized in that: The standard geographic spatial information adopts a topographic map with a scale of not less than 1:10,000 or an orthophoto with a spatial resolution of not less than 1 meter.

6. The method for improving geometric accuracy of historical geography maps according to claim 1, characterized in that: The single-band intensity image is obtained by weighted averaging of the multi-band images within the spectral range, or by using the first principal component of the multi-band images within the spectral range, or by using the first component of the Gram-Schmidt transform of the multi-band images within the spectral range.

7. The method for improving geometric accuracy of historical geography maps according to claim 1, characterized in that: When performing histogram matching, the single-band intensity image needs to be upsampled to the same resolution as the geometrically corrected high-resolution panchromatic image, and then the two need to be histogram matched so that the cumulative histogram of the high-resolution panchromatic image approaches the upsampled single-band intensity image. Finally, the two types of images after histogram matching are fused using the Gram-Schmidt algorithm to generate a historical surface benchmark color image.

8. The method for improving geometric accuracy of historical geography maps according to claim 2, characterized in that: When geometrically correcting and correcting the patch content of the historical geological map based on the two-dimensional reference base map and the three-dimensional reference base map, it is necessary to select the same-name features of the same period as the control points on the reference base map and the historical geological map, perform geometric correction on the historical geological map, and then compare with the reference base map to correct and improve the soil patch boundaries and attribute content on the geometrically corrected historical geological map. Finally, the final map is output after topological relationship checking and correction.

9. A computer-readable storage medium, characterized in that The storage medium stores a computer program, and when the computer program is executed by the processor, the method for improving the geometric accuracy correction of historical geographical maps as described in any one of claims 1 to 7 is implemented.

10. A computer electronic device, characterized in that: including memory and processor; The memory is used to store computer programs; The processor is used to implement the method for improving the geometric accuracy correction of historical geographical maps as described in any one of claims 1 to 7 when executing the computer program.

Images