Remote sensing image rapid selection method for natural resource normalized monitoring

Abstract

The invention relates to the field of remote sensing image processing, and particularly discloses a natural resource normalized monitoring-oriented remote sensing image rapid selection method, which comprises the following steps of S1, acquiring multi-source remote sensing image data, and performing preprocessing to obtain preprocessed image data containing semantic annotation; s2, carrying out coordinate conversion on the preprocessed image data containing the semantic annotation, obtaining feature parameters of the preprocessed image data containing the semantic annotation by adopting a coordinate mapping model based on deep learning, inputting the feature parameters into the coordinate mapping model for targeted correction, inputting an image data path, and outputting the image data path; reading a field name of a data attribute table in the image data path; s3, for the image data after coordinate conversion, firstly inputting data paths of a front batch of falling images and a rear batch of falling images, and reading the front batch of falling image data and the rear batch of falling image data; by adopting the technical scheme of the invention, the image data quality can be improved, the multi-source image coordinate mapping error is reduced, and the result reliability is ensured.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing technology, and in particular to a method for rapid selection of remote sensing images for routine monitoring of natural resources. Background Technology

[0002] In the field of routine natural resource monitoring, the accurate selection and efficient processing of remote sensing images are crucial for extracting change patches and supporting precise management of natural resources. Traditional monitoring methods rely on manpower, manually extracting change patches one by one. This approach has limitations, being done only once a year, resulting in weak business collaboration, low efficiency, and an inability to meet the real-time needs of farmland protection, early detection, early warning, early prevention of illegal land use, and ecological restoration.

[0003] Based on the above needs, Chinese Patent Publication No. CN118012863A discloses a method for rapid selection of remote sensing images for routine monitoring of natural resources. By acquiring remote sensing image data, coordinate transformation, image selection, and extraction of change patches, it improves the efficiency problem of traditional manual methods. However, it still has key defects due to the limitations of the technical design concept.

[0004] During use, due to the lack of preprocessing of remote sensing image data and the failure to consider the potential noise and radiometric imbalance of multi-source images (such as commercial and public satellite images), the default input image quality meets the standards, which makes subsequent coordinate transformation and image selection susceptible to interference from data quality.

[0005] Furthermore, the coordinate transformation uses a conventional process (reading the attribute table and generating a txt file). Because it lacks an optimized model designed for the coordinate mapping accuracy requirements of large-scale, multi-source imagery, and only implements the transformation using basic operations, it contains significant errors. In addition, the above scheme only addresses threshold setting and boundary correction for changing features extraction, without establishing an accuracy verification mechanism or designing result feedback iteration logic. This makes it impossible to confirm whether the extraction results meet the standards; if the accuracy is insufficient, manual backtracking is required, which increases costs.

[0006] Therefore, there is an urgent need for a rapid remote sensing image selection method for routine monitoring of natural resources that can improve image data quality, reduce coordinate mapping errors of multi-source images, and ensure the reliability of results. Summary of the Invention

[0007] This invention provides a method for rapid selection of remote sensing images for routine monitoring of natural resources, which can improve image data quality, reduce coordinate mapping errors of multi-source images, and ensure the reliability of results.

[0008] To solve the above-mentioned technical problems, this application provides the following technical solution: A method for rapid selection of remote sensing images for routine monitoring of natural resources includes the following: S1 acquires multi-source remote sensing image data, performs preprocessing, and obtains preprocessed image data with semantic annotations; S2 performs coordinate transformation on the preprocessed image data with semantic annotations. It adopts a deep learning-based coordinate mapping model to obtain the feature parameters of the preprocessed image data with semantic annotations, and inputs the feature parameters into the coordinate mapping model for targeted correction. Then, it inputs the image data path, reads the field names of the data attribute table in the image data path, and then creates a new blank file and writes the image result data after coordinate transformation corrected by the coordinate mapping model into the blank file to obtain the image data after coordinate transformation. S3 first inputs the data paths of the previous and subsequent batches of image data after coordinate transformation, and reads the previous and subsequent batches of image data; extracts ground feature points from the previous and subsequent batches of image data; matches the extracted ground feature points using the SIFT feature matching algorithm, and selects ground feature points with a matching degree higher than a preset threshold as reference control points; determines the operation range by combining the reference control points with data production requirements; based on the operation range, the image composite map is divided into blocks according to the overlapping areas of the preceding and following temporal images to obtain the image overlapping areas; after cleaning, the cleaned image composite map field is obtained; and a weighted fusion algorithm is used to process the image composite map to obtain the effective monitoring image composite map. S4 uses an improved U-Net semantic segmentation model to extract change patches from the effective monitoring images and outputs a dataset of key monitoring change patches that meets the accuracy requirements.

[0009] The basic principle and beneficial effects of this method are as follows: This method solves the core pain points of remote sensing image processing in the routine monitoring of natural resources through a four-stage progressive process of data preprocessing, accurate coordinate transformation, high-quality image selection, and reliable change patch extraction.

[0010] The S1 multi-source image preprocessing and semantic annotation principle addresses potential noise interference and radiometric imbalances in multi-source remote sensing images (such as commercial and public satellite images). By preprocessing, it eliminates data quality defects, laying a high-quality data foundation for subsequent stages. Simultaneously, semantic annotation is introduced to distinguish between natural changes, permissible restoration changes, and key monitored abnormal changes based on differences in ground features. This enables preliminary screening of monitoring targets and avoids non-key change data occupying subsequent processing resources.

[0011] S2 deep learning coordinate transformation breaks through the limitations of traditional one-size-fits-all coordinate transformation. By extracting feature parameters (such as land cover type and image clarity) from semantically labeled images, it makes targeted corrections to the deep learning coordinate mapping model, enabling the model to adapt to images with different data characteristics, establish more accurate coordinate mapping relationships, and reduce coordinate deviations caused by differences in data characteristics of multi-source images.

[0012] In the S3 image selection and fusion process, based on the previous and subsequent batches of map data, the SIFT feature matching algorithm is used to capture high-confidence ground feature points (such as road intersections and building corners) as reference control points. Combined with production requirements such as monitoring area boundaries and land type accuracy, the operation scope is locked to ensure that image selection focuses on the core monitoring area. Then, through block segmentation, cleaning, and weighted fusion, redundant and overlapping images are removed and the block stitching effect is optimized to form a complete and high-quality effective monitoring image combination map.

[0013] In the S4 improved U-Net change patch extraction, the improved U-Net semantic segmentation model is used to process the effective monitoring image combined with the map. This model has a stronger ability to capture ground feature details and can accurately identify the change areas in the images before and after the time phase. At the same time, through the built-in accuracy verification logic, it is ensured that the output dataset of key monitoring change patches meets the monitoring accuracy requirements and avoids invalid data output.

[0014] The S1 preprocessing stage eliminates noise interference from multi-source images (such as salt-and-pepper noise and Gaussian noise) and balances the radiometric differences between different satellite images (such as brightness deviations caused by different sensor performance between commercial and public satellites), thereby improving the signal-to-noise ratio of the image data. Semantic annotation further filters out core monitoring data, reduces the interference of non-key change data on subsequent processing, and improves the utilization rate of effective data.

[0015] The deep learning coordinate mapping model in S2, which is based on image feature parameter correction, can adapt to images of different land cover types (such as farmland and building land) and different resolutions. Compared with traditional coordinate transformation methods, the coordinate mapping error is controlled within 0.3-0.5 pixels, reducing the error and effectively solving the coordinate misalignment problem when stitching multi-source images.

[0016] S3 uses benchmark control points to lock the work area and cleans up redundant images to ensure the integrity and relevance of the combined imagery map. S4 improves the U-Net model to identify change areas with higher accuracy. Combined with the accuracy verification mechanism, it stabilizes the accuracy of extracting key monitoring change patches and avoids misjudgments and omissions caused by poor data quality and coordinate deviations, providing reliable data support for subsequent on-site verification.

[0017] On the one hand, S1 semantic annotation reduces the amount of non-critical data processing, shortening the processing time for subsequent coordinate transformation and image selection; on the other hand, S3 automated feature matching and block fusion processes replace traditional manual screening, avoiding efficiency bottlenecks caused by manpower-intensive methods, shortening the entire process cycle from image acquisition to change patch output, and meeting the needs of routine monitoring of natural resources.

[0018] For different monitoring scenarios (such as plain farmland, urban built-up areas, and mountainous forest land), the S2 coordinate mapping model can be adapted to the characteristics of the land cover by modifying the feature parameters (such as building land with clear boundaries, the model can increase the coordinate correction strength; forest land with complex textures, the model can optimize the detail capture capability). The S3 weighted fusion algorithm can adjust the block weight according to the land cover type (such as building land with higher splicing weight to ensure continuous boundaries), so that the method can stably output high-quality results in diverse geographical environments, breaking the limitation of traditional methods to adapt to a single scenario.

[0019] The output dataset of key monitoring change patches is not only highly accurate, but also associates change types (such as illegal land use and ecological restoration areas) through S1 semantic annotation. It can be directly connected to natural resource supervision business systems (such as farmland protection early warning platform and illegal land use verification system) to achieve rapid transformation from data to decision-making, accurately locate problem areas, and improve the level of intelligence in natural resource management.

[0020] In summary, this invention achieves the goals of improving image data quality, reducing coordinate mapping errors in multi-source images, and ensuring the reliability of results.

[0021] Furthermore, the preprocessing of S1 includes sequentially performing denoising processing to eliminate image noise interference, radiometric correction processing to balance image radiometric information, and then using a semantic understanding model to perform semantic annotation on the preprocessed image data to distinguish between natural changes caused by time / environment, changes in natural resources that can be restored, and abnormal changes that need to be monitored.

[0022] Further, in S3, the ground feature points include road intersections, building corners, and topographic landmarks. The extracted ground feature points are matched using the SIFT feature matching algorithm, and those with a matching degree higher than a preset threshold are selected as baseline control points. The operational scope is determined by combining the baseline control points with data production requirements, including monitoring area boundaries and land type accuracy requirements. Based on the operational scope, the image composite map is divided into blocks according to the overlapping areas of preceding and following temporal images to obtain overlapping image areas. The overlapping image areas are cleaned according to selection rules: within the same batch of images, images with a temporal difference greater than a first preset duration are selected in order of temporal phase, resolution, and cloud cover degree; images with a temporal difference within a first preset duration are selected in the same order; images from different batches are selected with a temporal interval greater than a second preset duration. After cleaning, a cleaned image composite map field is obtained. A weighted fusion algorithm is used to fuse adjacent blocks corresponding to the cleaned image composite map field to ensure the continuity of ground feature features at the block splicing points, resulting in an effective monitoring image composite map.

[0023] Furthermore, in S4, combining the semantic annotations obtained in S1, the initial change patches requiring key monitoring are extracted based on the effective production range of the monitoring images and the images from the preceding and following time periods. A grayscale difference threshold and a shape constraint threshold are set, and the initial change patches are post-processed. This post-processing includes removing small patches with areas smaller than the shape constraint threshold and merging adjacent patches with the same land use type. A polynomial smoothing algorithm is used to correct the boundaries of the post-processed initial change patches, reducing jagged lines, resulting in a change patch dataset. A confusion matrix is ​​used to verify the accuracy of the change patch dataset, calculating the overall accuracy and Kappa coefficient. If the overall accuracy is lower than a first preset value and the Kappa coefficient is lower than a second preset value, the process returns to S3 for image selection and processing again. If the overall accuracy is not lower than the first preset value and the Kappa coefficient is not lower than the second preset value, a key monitoring change patch dataset that meets the accuracy requirements is output.

[0024] Furthermore, the semantic annotation process of the semantic understanding model in S1 employs a multi-feature fusion semantic classification algorithm, including: S1-1 Extracts three types of core features from the preprocessed image data, including ground cover spectral features. Texture features Temporal variation characteristics Among them, the spectral characteristics of ground features The reflectance of each pixel in the image is obtained by calculating the average reflectance in the four bands: red (R), green (G), blue (B), and near-infrared (NIR). , , , , These represent the pixel reflectance values ​​for the corresponding wavelength bands, ranging from 0 to 1; texture features. Contrast calculated using the gray-level co-occurrence matrix With entropy The weighted sum, i.e. , Indicates the degree of lightness and darkness difference in the texture of ground features, with a value range of 0-255. Indicates the degree of disorder in texture information, with a value ranging from 0 to 8; temporal variation characteristics. The mean value of the absolute values ​​of the pixel grayscale differences between the current image and reference images of the same area within the previous three months is obtained, i.e. , This indicates the number of reference images used in the calculation. , Indicates the current image number The grayscale value of the pixel corresponding to each reference time period. Indicates the first The grayscale value of a pixel corresponding to a reference image; S1-2 inputs the three core features into the classification layer of the semantic understanding model, and uses a classification function. Output semantic annotation results ; in For feature vectors, , , These are the weight matrices from the model input layer to hidden layer 1, from hidden layer 1 to hidden layer 2, and from hidden layer 2 to the output layer, with dimensions of 64×3, 32×64, and 3×32, respectively. , , These are the bias vectors for the corresponding layers, with dimensions of 64×1, 32×1, and 3×1, respectively. It is the Sigmoid activation function. A 3D vector ( ), , , These represent the probabilities of the corresponding pixel belonging to natural changes caused by time / environment, changes that can be repaired by natural resources, and abnormal changes that require key monitoring, respectively. The values ​​range from 0 to 1, and the category with the highest probability is the final semantic annotation category.

[0025] Furthermore, when performing targeted corrections on the deep learning-based coordinate mapping model in S2, an adaptive correction algorithm for feature parameters is employed, including: S2-1 Extracts two types of key feature parameters from the semantically labeled preprocessed image data: the proportion of land cover types. With image sharpness parameters Among them, the proportion parameter of land cover type This is the area proportion vector of different land cover types in the image, i.e. , , , , These represent the area percentages of cultivated land, construction land, forest land, and water area, respectively, with values ​​ranging from 0 to 1. ; Image sharpness parameters It is obtained by calculating the mean value of the gradient values ​​at the image edges, that is , , These are the width and height of the image, respectively. , They respectively represent the images in Sobel gradient values ​​at a pixel along the horizontal and vertical directions; S2-2 Constructing the Correction Factor and The translation parameters of the coordinate mapping model are respectively With scaling parameters Make corrections, and then adjust the translation parameters. Corrected scaling parameters ,in This is the adjustment step size for the model's default translation parameters, with a value range of 0.1-0.5 pixels; The adjustment step size for the model's default scaling parameters, ranging from 0.001 to 0.005; correction factor. Based on the proportion parameter of land cover type Calculated, i.e. The higher the proportion of land used for construction, The larger the value, the clearer the building site boundaries, requiring a larger translation correction to match the coordinate reference; correction factor Image sharpness parameters Calculated, i.e. ,when hour, Follow Reduce and lower the resolution to avoid coordinate deviations caused by over-scaling of low-resolution images; when hour, The value is kept stable at 1.2 to ensure the accuracy of coordinate mapping for high-definition images.

[0026] This method can adapt to the coordinate correction needs of different monitoring areas (such as concentrated farmland areas and urban building areas) by using the land cover type proportion parameter, and can also avoid coordinate deviation caused by improper correction of low-quality images by using the image clarity parameter, so that the coordinate transformation error is controlled within 0.5 pixels. At the same time, it can adapt to the coordinate transformation needs of multi-source remote sensing images (commercial satellite and public welfare satellite images).

[0027] Furthermore, the extraction of ground feature points in S3 employs a multi-scale feature enhancement extraction algorithm, including: S3-1 constructs image pyramids at three scales for the previous batch of image data and the subsequent batch of image data: scale 1: original resolution, scale 2: 0.5 times the original resolution, and scale 3: 0.25 times the original resolution. Gaussian filtering is applied at each scale, and the kernel size increases with the scale: 3×3 for scale 1, 5×5 for scale 2, and 7×7 for scale 3. S3-2 uses the characteristic response function Calculate the feature response value for each pixel, where For pixel coordinates, The Gaussian standard deviation at the current scale (scale 1) Scale 2 Scale 3 ), scale factor ( ), The pixel grayscale value after Gaussian filtering at the corresponding scale; S3-3 Set the response threshold The value ranges from 10 to 20, and is adaptively adjusted based on image contrast. Response values ​​are then filtered. The pixels are used as candidate feature points, and non-local maxima in adjacent candidate points are eliminated by non-maximum suppression to finally obtain the ground feature points.

[0028] This method can cover features of different sizes through multi-scale pyramids, such as capturing road intersections at a small scale and capturing terrain landmarks at a large scale. It can also reduce the interference of noise on feature extraction through Gaussian filtering and non-maximum suppression, thereby improving the accuracy of feature point extraction and reducing the mismatch rate of subsequent SIFT matching.

[0029] Furthermore, the weighted fusion algorithm in S3 adopts an adaptive weighting strategy based on land cover type, including: S3-4 Combines the cleaned image with adjacent blocks corresponding to the map fields, and extracts the land cover type proportion vector for each block. , , representing two adjacent blocks respectively; -Proportion of arable land -Proportion of building land, - Forest land ratio, - Percentage of water area, with a value range of 0-1; S3-5 uses the weight calculation function Determine the fusion weight for each block. Among them, building land has the highest weight because the boundaries of building land require high-precision splicing, while water area has the lowest weight because water area has uniform grayness and high splicing tolerance. S3-6 For overlapping pixels in adjacent blocks, the fusion formula is used... Calculate the grayscale value after fusion, where , Block 1 and Block 2 are respectively in The grayscale value of a pixel.

[0030] This strategy can dynamically adjust the weights according to the stitching requirements of different land cover types. For example, the weight of blocks in areas with concentrated building land is higher to ensure continuous boundaries. It can also avoid obvious gray-scale jumps at the block stitching points through weighted averaging, thereby improving the stitching accuracy of the combined monitoring image map. At the same time, it can adapt to the image quality requirements of different monitoring scenarios such as farmland and building land.

[0031] Furthermore, the polynomial smoothing algorithm in S4 employs a boundary curvature constraint smoothing strategy, including: S4-1 extracts the boundary pixel sequence from the post-processed initial changed patch boundary. , Calculate the curvature of three adjacent pixels to represent the number of boundary pixels. , ,in Indicates the first The curvature of each boundary segment; the greater the curvature, the more severe the boundary bending. S4-2 sets the curvature threshold. The value ranges from 0.05 to 0.1, and it is adaptively adjusted according to the size of the patch. For the severely curved sections, a cubic polynomial was used for fitting: , ( , , (Solved using boundary smoothing constraints); For the smooth sections, keep the original boundary pixels unchanged; S4-3 stitches the fitted boundary pixels with the original smooth segment pixels to obtain the corrected patch boundary.

[0032] This strategy can preserve the key morphological features of the map patch boundary through curvature constraints, avoid excessive smoothing that leads to map patch distortion, effectively eliminate jagged lines, improve the smoothness of the boundary of the changing map patch, and improve the matching degree between the subsequent changing map patch and the field verification.

[0033] Furthermore, the accuracy verification in S4 employs a hierarchical accuracy weighted calculation method, including: S4-4 divides the changed patch dataset into three sub-datasets according to the semantic annotation in S1: natural change sub-dataset Repairing changed subsets of data Abnormal change subset Calculate the confusion matrix for each subset of data to obtain the overall accuracy within each class. , , The value range is 0-1; S4-5 sets the weight coefficients for the three types of subsets based on the weight requirements of the monitoring business. , , Abnormal changes have the highest weight and require close monitoring; natural changes have the lowest weight because misjudgment has a smaller impact, and are weighted using a weighted formula. Calculate the overall accuracy; S4-6 uses the same weighting coefficients to calculate the stratified Kappa coefficients. , These are the Kappa coefficients for the three subsets of data; S4-7if and If the accuracy meets the standard, then the process is considered to be satisfactory; otherwise, return to S3 for reprocessing.

[0034] This method can not only evaluate the extraction accuracy of different types of change patches in a hierarchical manner, avoiding the problem of low accuracy of key monitoring categories being masked by a single accuracy index, but also highlight the accuracy requirements of abnormal changes through business weights, making the accuracy control of key monitoring categories more precise, while reducing the interference of misjudgment of non-key changes on the overall accuracy assessment. Attached Figure Description

[0035] Figure 1 This is a flowchart illustrating an embodiment of a method for rapid selection of remote sensing images for routine monitoring of natural resources. Detailed Implementation

[0036] The following detailed description illustrates the specific implementation method: A method for rapid selection of remote sensing images for routine monitoring of natural resources (e.g.) Figure 1 As shown), it includes the following: S1 acquires multi-source remote sensing image data, performs preprocessing, and obtains preprocessed image data with semantic annotations; S2 performs coordinate transformation on the preprocessed image data with semantic annotations. It adopts a deep learning-based coordinate mapping model to obtain the feature parameters of the preprocessed image data with semantic annotations, and inputs the feature parameters into the coordinate mapping model for targeted correction. Then, it inputs the image data path, reads the field names of the data attribute table in the image data path, and then creates a new blank file and writes the image result data after coordinate transformation corrected by the coordinate mapping model into the blank file to obtain the image data after coordinate transformation. S3 first inputs the data paths of the previous and subsequent batches of image data after coordinate transformation, and reads the previous and subsequent batches of image data; extracts ground feature points from the previous and subsequent batches of image data; matches the extracted ground feature points using the SIFT feature matching algorithm, and selects ground feature points with a matching degree higher than a preset threshold as reference control points; determines the operation range by combining the reference control points with data production requirements; based on the operation range, the image composite map is divided into blocks according to the overlapping areas of the preceding and following temporal images to obtain the image overlapping areas; after cleaning, the cleaned image composite map field is obtained; and a weighted fusion algorithm is used to process the image composite map to obtain the effective monitoring image composite map. S4 uses an improved U-Net semantic segmentation model to extract change patches from the effective monitoring images and outputs a dataset of key monitoring change patches that meets the accuracy requirements.

[0037] The preprocessing of S1 includes sequentially performing denoising processing to eliminate image noise interference, radiometric correction processing to balance image radiometric information, and then using a semantic understanding model to perform semantic annotation on the preprocessed image data to distinguish between natural changes caused by time / environment, changes in natural resources that can be restored, and abnormal changes that need to be monitored.

[0038] In S3, the ground feature points include road intersections, building corners, and topographic landmarks. The extracted ground feature points are matched using the SIFT feature matching algorithm, and those with a matching degree higher than a preset threshold are selected as baseline control points. The operational scope is determined by combining the baseline control points with data production requirements, including monitoring area boundaries and land type accuracy requirements. Based on the operational scope, the image composite map is divided into blocks according to the overlapping areas of consecutive temporal images to obtain overlapping image areas. The overlapping image areas are cleaned according to selection rules: within the same batch of images, images with a temporal difference greater than a first preset duration are selected in order of temporal phase, resolution, and cloud cover degree; images with a temporal difference within the first preset duration are selected in the same order; images from different batches are selected with a temporal interval greater than a second preset duration. After cleaning, a cleaned image composite map field is obtained. A weighted fusion algorithm is used to fuse adjacent blocks corresponding to the cleaned image composite map field to ensure the continuity of ground feature features at the block splicing points, resulting in an effective monitoring image composite map.

[0039] In S4, combining the semantic annotations obtained in S1, the initial change patches requiring key monitoring are extracted based on the effective production range of the monitoring images and the images from the preceding and following time periods. A grayscale difference threshold and a shape constraint threshold are set, and the initial change patches are post-processed. This post-processing includes removing small patches with areas smaller than the shape constraint threshold and merging adjacent patches with the same land use type. A polynomial smoothing algorithm is used to correct the boundaries of the post-processed initial change patches, reducing jagged lines, resulting in a change patch dataset. A confusion matrix is ​​used to verify the accuracy of the change patch dataset, calculating the overall accuracy and Kappa coefficient. If the overall accuracy is lower than a first preset value and the Kappa coefficient is lower than a second preset value, the process returns to S3 for image selection and processing again. If the overall accuracy is not lower than the first preset value and the Kappa coefficient is not lower than the second preset value, a key monitoring change patch dataset that meets the accuracy requirements is output.

[0040] In practical use: In this embodiment, a certain area is designated as the monitoring area, and routine monitoring of natural resources is carried out in the first quarter to achieve dynamic monitoring of farmland protection, illegal land use supervision, and ecological restoration.

[0041] First, data preparation is carried out by acquiring multi-source remote sensing imagery data of the monitoring area, assuming it includes imagery data from a commercial satellite (resolution 0.8m, image taken on, for example, March 10, 2024) and a public satellite (resolution 2m, image taken on, for example, March 5, 2024). At the same time, the previous batch of mapping data (for example, the mapping data of a natural resource change survey of a certain area in December 2023, in Shapefile format), the subsequent batch of mapping data (orthophoto mapping data of the monitoring area in March 2024, in Shapefile format), and the boundary vector data of the monitoring area (the percentage of boundary delineation in the area in 2024) are acquired.

[0042] Then, the equipment and software were prepared, using a server configured with an Intel Xeon Gold 6338 processor, 128GB of memory, and 2TB of SSD storage; the software environment included ArcGIS Pro 3.1 (for image coordinate transformation and map data processing), TensorFlow 2.15 (for deep learning model deployment, including semantic understanding models and coordinate mapping models), and Python 3.9 (for algorithm programming implementation).

[0043] S1: Acquire multi-source remote sensing image data and preprocess it to obtain preprocessed image data with semantic annotations; Multi-source image preprocessing and denoising: For salt-and-pepper noise (noise pixels account for about 3%) in a commercial satellite image, a median filtering algorithm was used, with the filter window size set to 3×3; for Gaussian noise (standard deviation σ=1.2) in a public satellite image, a Gaussian filtering algorithm was used, with the filter kernel size set to 5×5. After denoising, the signal-to-noise ratio of the two types of images was improved from 28dB to 42dB. Radiometric correction was performed using ENVI 5.6 software. For a commercial satellite image, absolute radiometric correction based on the atmospheric radiative transfer model (6S model) was used. The atmospheric parameters at the time of image capture were input (assuming an aerosol optical thickness of 0.25 and an atmospheric water vapor content of 2.5 g / cm²), and the image grayscale values ​​were converted into surface reflectance. For a public satellite image, relative radiometric correction (histogram matching method) was used. Using the commercial satellite image as a reference, the mean deviation of the radiance between the two types of images was controlled within 5%.

[0044] Semantic annotation (using a multi-feature fusion semantic classification algorithm) Step S1-1: Extract 3 types of core features; Spectral characteristics of ground features Select any pixel in the image (e.g., coordinates (1200, 800)). Its red band (R) reflectance is 0.32, green band (G) is 0.28, blue band (B) is 0.25, and near-infrared band (NIR) is 0.65. According to the formula... Assuming the result is .

[0045] Texture features Calculate the gray-level co-occurrence matrix (distance 1, angle 0°) of the 3×3 neighborhood of the pixel to obtain the contrast ratio. Entropy According to the formula Assuming the result is .

[0046] Temporal variation characteristics Reference images of the same area were selected from the previous three months (January and February 2024 images from a public satellite). Supplementing images from December 2023 This pixel has the following grayscale value in the current image. (January reference) (February reference) (December reference), corresponding to the grayscale value of the reference image. , , According to the formula Assuming the result is .

[0047] Step S1-2: Semantic classification and annotation, converting feature vectors Input semantic understanding model, model weight matrix (64×3) (32×64) (3×32) Pre-training optimization is performed using stochastic gradient descent, with bias vectors... (64×1) (32×1) (3×1) The initial value is set to 0.1; after calculation by the Sigmoid activation function, the semantic annotation result is output. ,in (The abnormal changes that need to be monitored are the most likely), so this pixel is marked as an abnormal change that needs to be monitored (subsequent verification confirmed that the area was a building site built without prior approval).

[0048] The final result is preprocessed image data with semantic annotations (in TIFF format). The annotation categories include natural changes caused by time / environment (such as areas with seasonal vegetation withering, accounting for about 15%), changes that can be restored by natural resources (such as mine ecological restoration areas, accounting for about 8%), and abnormal changes that need to be monitored (such as areas with illegal land use and illegal construction, accounting for about 5%).

[0049] S2: Perform coordinate transformation on the preprocessed image data with semantic annotations to obtain the image data after coordinate transformation; The feature parameter adaptive correction algorithm is adopted, and the specific steps are as follows: Step S2-1: Extract key feature parameters; Ground feature type proportion parameter For the preprocessed image data with semantic annotations, the monitored area includes 120 km² of cultivated land, 80 km² of built-up land, 60 km² of forest land, and 40 km² of water area, totaling 300 km². According to the formula... Assuming the result is , , , ,Right now .

[0050] Image sharpness parameters Width of a certain commercial satellite image Pixels, height Pixel, for any pixel in the image ( ), calculate its horizontal Sobel gradient value Vertical direction According to the formula The sum of the gradient values ​​of the entire image is calculated as follows: Assuming the result is (Due to the high resolution of a certain commercial satellite image, the actual sharpness parameters are...) (This is an example calculation).

[0051] Step S2-2: Coordinate mapping model correction and transformation; Construct the correction factor according to the formula Assuming the result is Image clarity According to the formula Assuming the result is .

[0052] Correcting model parameters: Default translation parameters for the model Pixels, Adjust Step Size Pixels, corrected translation parameters Pixels; Default scaling parameters Adjust step size Assuming the scaling parameters are corrected .

[0053] The coordinate transformation was performed by inputting the path of the preprocessed image data with semantic annotations (D:\Data\Preprocessed\) in ArcGIS Pro, reading the field names of the data attribute table (name, scenetime, cloudperce, imagegsd), creating a new txt file (D:\Data\ImageNames.txt) and writing the image names ("20240310.tif" "GF6_20240305.tif"); inputting the coordinate-transformed image after model correction, reading the txt file to determine the actual map area of ​​the image (the map area of ​​a commercial satellite image: longitude XX-XX, latitude XX-XX), and outputting the image data converted to the geodetic coordinate system (in TIFF format). After inspection, the coordinate transformation error was controlled within 0.4 pixels, which is lower than the error of traditional transformation methods.

[0054] S3: Process the image data after coordinate transformation to obtain an effective monitoring image combination map; extract ground feature points using a multi-scale feature enhancement extraction algorithm; Step S3-1: Construct image pyramids and apply Gaussian filtering. For the previous batch of map data (December 2023) and the subsequent batch of map data (March 2024), construct image pyramids at three scales: Scale 1 (original resolution 1m), Scale 2 (0.5m), and Scale 3 (0.25m). Apply Gaussian filtering at each scale, with a 3×3 filter kernel for Scale 1. ), Scale 2 cores 5×5 ( ), 3 cores 7×7 ( After filtering, the image noise is reduced.

[0055] Step S3-2: Calculate the feature response value and take the pixels of road intersections from the subsequent batch of map data. ), scale 1 , Assuming the pixel grayscale value after Gaussian filtering , , According to the formula Assuming the result is ; Scale 2 Assuming the result is ; Scale 3 Assuming the result is .

[0056] Step S3-3: Filter ground feature points and set a response threshold based on image contrast (contrast of 80 for a certain commercial satellite image). Filter response values The pixels are used as candidate feature points (such as building corner points). ); By using 3×3 neighborhood nonmaxima suppression, nonlocal maxima (such as adjacent pixels) are eliminated. (Points less than 18 were removed), and a total of 2,000 feature points were finally extracted, including 800 road intersections, 700 building corner points, and 500 terrain landmarks (such as mountain tops).

[0057] The baseline control points were determined and the work area was delineated. The SIFT feature matching algorithm was used to match the feature points of the land cover data in the previous and subsequent batches. A matching degree threshold of 0.85 was set (using Euclidean distance as a measure, a distance less than 30 is considered a successful match). 1600 feature points with a matching degree higher than 0.85 were selected as baseline control points.

[0058] Based on the baseline control points and data production requirements, the operation area was determined to be XX-XX″ east longitude and XX-XX north latitude, with an area of ​​approximately 280 km² (excluding edge areas without data) through buffer analysis (buffer distance 50m) and overlay analysis in ArcGIS.

[0059] Block processing: Based on the scope of the operation, the combined image (a mosaic image of a commercial satellite and a public welfare satellite) is divided into blocks according to the overlapping area of ​​the preceding and following time phase images (overlap rate 30%). Each block is 1000×1000 pixels in size, resulting in a total of 300 blocks, of which 120 are overlapping area blocks.

[0060] The image cleaning process is performed with two preset durations: 3 days for the first preset duration and 15 days for the second. Within the same batch (e.g., images from a commercial satellite on March 10th and 12th, 2024), if the time phase difference is 2 days (≤3 days), the image is selected based on resolution, time phase, and cloud cover: the March 10th image has a resolution of 0.8m (higher than the March 12th image's 1m) and a cloud cover of 5% (lower than the March 12th image's 8%), therefore the March 10th image is retained. In different batches (e.g., images from a commercial satellite on March 10th and a public interest satellite on March 5th), if the time phase difference is 5 days (≤15 days), the public interest satellite image on March 5th is not selected; instead, the commercial satellite image on March 10th is used. After cleaning, 30 overlapping areas that do not meet the requirements are removed, resulting in cleaned image-combined map fields (including block numbers, land cover type percentages, etc.).

[0061] Weighted fusion employs an adaptive weighting strategy based on land cover type. Step S3-4: Extract the land cover type proportion vector, and take adjacent block 1 (numbered A01) and block 2 (numbered A02). The proportion of cultivated land in block 1 is... Construction land ratio forest land ratio water area ratio ,Right now ; Proportion of cultivated land in block 2 Construction land ratio forest land ratio water area ratio ,Right now .

[0062] Step S3-5: Calculate the fusion weights according to the formula. Assuming the result is block 1 weight Block 2 weight (Because Block 2 has a higher proportion of building land, it has a greater weight).

[0063] Steps S3-6, pixel fusion, for overlapping pixels ( ), grayscale value of block 1 2. Grayscale value of block 2 According to the formula Assuming the result is the merged grayscale value After fusion, the grayscale change at the block splicing point decreased from 15 to 3, and the ground features remained continuous without any faults.

[0064] The final result is an effective monitoring image combination map. The boundaries of ground features in the combination map are clear, with no overlapping or redundant areas, which meets the requirements for subsequent change patch extraction.

[0065] S4: Extract change patches and output a dataset of key monitored change patches. Initial change patch extraction: Load the improved U-Net semantic segmentation model (3 input layer channels, 3 output layer channels, corresponding to 3 change types) into TensorFlow. Input the effective monitoring image combined with images from two periods (December 2023 and March 2024). Combine the semantic annotations from S1 (focusing on abnormal changes that need to be monitored) and set the change detection threshold (grayscale difference ≥ 20). Extract a total of 120 initial change patches, including 45 patches of farmland converted to building land, 30 patches of forest land converted to construction land, and 45 other types of change patches.

[0066] Post-processing and boundary correction of image features employ a boundary curvature constraint smoothing strategy. Step S4-1: Calculate the boundary curvature. Take a map patch of farmland converted to construction land (area 1500m²) and extract the boundary pixel sequence. ( ), for the 10th boundary segment ( ), , , , , , According to the formula Assuming the result is .

[0067] Step S4-2: Polynomial fitting and smoothing, setting curvature threshold. ,because A cubic polynomial was used to fit this segment: , ( ), through boundary smoothing constraints (fitting curve over and (and the first derivative is continuous), solving yields , , , , , ;right The gentle section (such as) , ), keeping the original boundary pixels unchanged.

[0068] Step S4-3, boundary stitching: stitch the fitted boundary pixels with the smooth segment pixels to obtain the corrected patch boundary. The jagged lines are reduced from 3 times per 10 pixels to 0.5 times, and the boundary smoothness is improved.

[0069] Post-processing of the patches involves setting a shape constraint threshold (minimum patch area 50m²) to remove 15 tiny patches with an area less than 50m²; merging adjacent patches with the same land use type (e.g., 3 adjacent farmland-to-construction patches are merged into 1 patch with an area of ​​800m²) to finally obtain a dataset of changed patches (containing 100 patches in GeoJSON format).

[0070] Accuracy verification employs a hierarchical accuracy weighted calculation method; Step S4-4: Divide the dataset into sub-datasets and calculate intra-class precision. Divide the changed patch dataset into three sub-datasets according to the semantic annotation in S1: natural change sub-dataset. (20 patches), Repairing the changed subset of the dataset (30 map patches), anomalous change subset (50 map features); Through on-site verification (using RTK positioning, the number of verification samples is 30% of the number of map features), the result is assumed to be the overall accuracy within the class. , , Kappa coefficient , , .

[0071] Step S4-5: Calculate the overall accuracy and set the weighting coefficients. , , According to the formula Assuming the result is .

[0072] Step S4-6: Calculate the stratified Kappa coefficients according to the formula. Assuming the result is .

[0073] Step S4-7: Accuracy Determination and Result Output: Because and If the accuracy is deemed satisfactory, a dataset of key monitoring change patches (containing 50 abnormal change patches, along with attribute information such as patch location, area, and change type) that meets the accuracy requirements will be output. This dataset can be directly output to the natural resources supervision platform for the investigation and handling of illegal land use and the supervision of ecological restoration.

[0074] This embodiment improves the signal-to-noise ratio and keeps radiometric differences within a low range after preprocessing multi-source images; semantic annotation accuracy is high, effectively filtering core monitoring data and reducing interference from non-key data. Coordinate mapping errors are small, meeting the accuracy requirements for multi-source image stitching and land cover identification. The efficiency of obtaining change patches from images is high. The accuracy of extracting key monitoring change patches is high.

[0075] The above are merely embodiments of the present invention. The invention is not limited to the fields covered by these embodiments. Commonly known structures and characteristics in the solutions are not described in detail here. Those skilled in the art are aware of all common technical knowledge in the field prior to the application date or priority date, are able to access all existing technologies in that field, and have the ability to apply conventional experimental methods prior to that date. Those skilled in the art can, under the guidance of this application, improve and implement this solution in combination with their own capabilities. Some typical known structures or methods should not be obstacles for those skilled in the art to implement this application. It should be noted that those skilled in the art can make several modifications and improvements without departing from the structure of the present invention. These should also be considered within the scope of protection of the present invention, and will not affect the effectiveness of the implementation of the present invention or the practicality of the patent. The scope of protection claimed in this application should be determined by the content of its claims, and the specific embodiments described in the specification can be used to interpret the content of the claims.

Claims

1. A method for rapid selection of remote sensing images for routine monitoring of natural resources, characterized in that, Includes the following: S1 acquires multi-source remote sensing image data, performs preprocessing, and obtains preprocessed image data with semantic annotations; S2 performs coordinate transformation on the preprocessed image data with semantic annotations. It adopts a deep learning-based coordinate mapping model to obtain the feature parameters of the preprocessed image data with semantic annotations, and inputs the feature parameters into the coordinate mapping model for targeted correction. Then, it inputs the image data path, reads the field names of the data attribute table in the image data path, and then creates a new blank file and writes the image result data after coordinate transformation corrected by the coordinate mapping model into the blank file to obtain the image data after coordinate transformation. S3 first inputs the data paths of the previous and subsequent batches of image data after coordinate transformation, and reads the previous and subsequent batches of image data; extracts ground feature points from the previous and subsequent batches of image data; matches the extracted ground feature points using the SIFT feature matching algorithm, and selects ground feature points with a matching degree higher than a preset threshold as reference control points; determines the operation range by combining the reference control points with data production requirements; based on the operation range, the image composite map is divided into blocks according to the overlapping areas of the preceding and following temporal images to obtain the image overlapping areas; after cleaning, the cleaned image composite map field is obtained; and a weighted fusion algorithm is used to process the image composite map to obtain the effective monitoring image composite map. S4 uses an improved U-Net semantic segmentation model to extract change patches from the effective monitoring images and outputs a dataset of key monitoring change patches that meets the accuracy requirements.

2. The method for rapid selection of remote sensing images for routine monitoring of natural resources according to claim 1, characterized in that, The preprocessing of S1 includes sequentially performing denoising processing to eliminate image noise interference, radiometric correction processing to balance image radiometric information, and then using a semantic understanding model to perform semantic annotation on the preprocessed image data to distinguish between natural changes caused by time / environment, changes in natural resources that can be restored, and abnormal changes that need to be monitored.

3. The method for rapid selection of remote sensing images for routine monitoring of natural resources according to claim 2, characterized in that, In S3, the ground feature points include road intersections, building corners, and topographic landmarks. The extracted ground feature points are matched using the SIFT feature matching algorithm, and those with a matching degree higher than a preset threshold are selected as baseline control points. The operational scope is determined by combining the baseline control points with data production requirements, including monitoring area boundaries and land type accuracy requirements. Based on the operational scope, the image composite map is divided into blocks according to the overlapping areas of consecutive temporal images to obtain overlapping image areas. The overlapping image areas are cleaned according to selection rules: within the same batch of images, images with a temporal difference greater than a first preset duration are selected in order of temporal phase, resolution, and cloud cover degree; images with a temporal difference within the first preset duration are selected in the same order; images from different batches are selected with a temporal interval greater than a second preset duration. After cleaning, a cleaned image composite map field is obtained. A weighted fusion algorithm is used to fuse adjacent blocks corresponding to the cleaned image composite map field to ensure the continuity of ground feature features at the block splicing points, resulting in an effective monitoring image composite map.

4. The method for rapid selection of remote sensing images for routine monitoring of natural resources according to claim 3, characterized in that, In S4, combining the semantic annotations obtained in S1, the initial change patches requiring key monitoring are extracted based on the effective production range of the monitoring images and the images from the preceding and following time periods. A grayscale difference threshold and a shape constraint threshold are set, and the initial change patches are post-processed. This post-processing includes removing small patches with areas smaller than the shape constraint threshold and merging adjacent patches with the same land use type. A polynomial smoothing algorithm is used to correct the boundaries of the post-processed initial change patches, reducing jagged lines, resulting in a change patch dataset. A confusion matrix is ​​used to verify the accuracy of the change patch dataset, calculating the overall accuracy and Kappa coefficient. If the overall accuracy is lower than a first preset value and the Kappa coefficient is lower than a second preset value, the process returns to S3 for image selection and processing again. If the overall accuracy is not lower than the first preset value and the Kappa coefficient is not lower than the second preset value, a key monitoring change patch dataset that meets the accuracy requirements is output.

5. The method for rapid selection of remote sensing images for routine monitoring of natural resources according to claim 4, characterized in that, The semantic annotation process of the semantic understanding model in S1 employs a multi-feature fusion semantic classification algorithm, including: S1-1 Extracts three types of core features from the preprocessed image data, including ground cover spectral features. Texture features Temporal variation characteristics Among them, the spectral characteristics of ground features The reflectance of each pixel in the image is obtained by calculating the average reflectance in the four bands: red (R), green (G), blue (B), and near-infrared (NIR). , , , , These represent the pixel reflectance values ​​for the corresponding wavelength bands, ranging from 0 to 1; texture features. Contrast calculated using the gray-level co-occurrence matrix With entropy The weighted sum, i.e. , Indicates the degree of lightness and darkness difference in the texture of ground features, with a value range of 0-255. Indicates the degree of disorder in texture information, with a value ranging from 0 to 8; temporal variation characteristics. The mean value of the absolute values ​​of the pixel grayscale differences between the current image and reference images of the same area within the previous three months is obtained, i.e. , This indicates the number of reference images used in the calculation. , Indicates the current image number The grayscale value of the pixel corresponding to each reference time period. Indicates the first The grayscale value of a pixel corresponding to a reference image; S1-2 inputs the three core features into the classification layer of the semantic understanding model, and uses a classification function. Output semantic annotation results ; in For feature vectors, , , These are the weight matrices from the model input layer to hidden layer 1, from hidden layer 1 to hidden layer 2, and from hidden layer 2 to the output layer, with dimensions of 64×3, 32×64, and 3×32, respectively. , , These are the bias vectors for the corresponding layers, with dimensions of 64×1, 32×1, and 3×1, respectively. It is the Sigmoid activation function. A 3D vector ( ), , , These represent the probabilities of the corresponding pixel belonging to natural changes caused by time / environment, changes that can be repaired by natural resources, and abnormal changes that require key monitoring, respectively. The values ​​range from 0 to 1, and the category with the highest probability is the final semantic annotation category.

6. The method for rapid selection of remote sensing images for routine monitoring of natural resources according to claim 5, characterized in that, When performing targeted corrections on the deep learning-based coordinate mapping model in S2, an adaptive correction algorithm for feature parameters is employed, including: S2-1 Extracts two types of key feature parameters from the semantically labeled preprocessed image data: the proportion of land cover types. With image sharpness parameters Among them, the proportion parameter of land cover type This is the area proportion vector of different land cover types in the image, i.e. , , , , These represent the area percentages of cultivated land, construction land, forest land, and water area, respectively, with values ​​ranging from 0 to 1. ; Image sharpness parameters It is obtained by calculating the mean value of the gradient values ​​at the image edges, that is , , These are the width and height of the image, respectively. , They respectively represent the images in Sobel gradient values ​​at a pixel along the horizontal and vertical directions; S2-2 Constructing the Correction Factor and The translation parameters of the coordinate mapping model are respectively With scaling parameters Make corrections, and then adjust the translation parameters. Corrected scaling parameters ,in This is the adjustment step size for the model's default translation parameters, with a value range of 0.1-0.5 pixels; The adjustment step size for the model's default scaling parameters, ranging from 0.001 to 0.005; correction factor. Based on the proportion parameter of land cover type Calculated, i.e. The higher the proportion of land used for construction, The larger the value, the clearer the building site boundaries, requiring a larger translation correction to match the coordinate reference; correction factor Image sharpness parameters Calculated, i.e. ,when hour, Follow Reduce and lower the resolution to avoid coordinate deviations caused by over-scaling of low-resolution images; when hour, The value is kept stable at 1.2 to ensure the accuracy of coordinate mapping for high-definition images.

7. The method for rapid selection of remote sensing images for routine monitoring of natural resources according to claim 6, characterized in that, The S3 method employs a multi-scale feature enhancement extraction algorithm when extracting ground feature points, including: S3-1 constructs image pyramids at three scales for the previous batch of image data and the subsequent batch of image data: scale 1: original resolution, scale 2: 0.5 times the original resolution, and scale 3: 0.25 times the original resolution. Gaussian filtering is applied at each scale, and the kernel size increases with the scale: 3×3 for scale 1, 5×5 for scale 2, and 7×7 for scale 3. S3-2 uses the characteristic response function Calculate the feature response value for each pixel, where For pixel coordinates, The Gaussian standard deviation at the current scale (scale 1) Scale 2 Scale 3 ), scale factor ( ), The pixel grayscale value after Gaussian filtering at the corresponding scale; S3-3 Set the response threshold The value ranges from 10 to 20, and is adaptively adjusted based on image contrast. Response values ​​are then filtered. The pixels are used as candidate feature points, and non-local maxima in adjacent candidate points are eliminated by non-maximum suppression to finally obtain the ground feature points.

8. A method for rapid selection of remote sensing images for routine monitoring of natural resources according to claim 7, characterized in that, The weighted fusion algorithm in S3 adopts an adaptive weighting strategy based on land cover type, including: S3-4 Combines the cleaned image with adjacent blocks corresponding to the map fields, and extracts the land cover type proportion vector for each block. , , representing two adjacent blocks respectively; -Proportion of arable land -Proportion of building land, - Forest land ratio, - Percentage of water area, with a value range of 0-1; S3-5 uses the weight calculation function Determine the fusion weight for each block. Among them, building land has the highest weight because the boundaries of building land require high-precision splicing, while water area has the lowest weight because water area has uniform grayness and high splicing tolerance. S3-6 For overlapping pixels in adjacent blocks, the fusion formula is used... Calculate the grayscale value after fusion, where , Block 1 and Block 2 are respectively in The grayscale value of a pixel.

9. A method for rapid selection of remote sensing images for routine monitoring of natural resources according to claim 8, characterized in that, The polynomial smoothing algorithm in S4 employs a boundary curvature constraint smoothing strategy, including: S4-1 extracts the boundary pixel sequence from the post-processed initial changed patch boundary. , Calculate the curvature of three adjacent pixels to represent the number of boundary pixels. , ,in Indicates the first The curvature of each boundary segment; the greater the curvature, the more severe the boundary bending. S4-2 sets the curvature threshold. The value ranges from 0.05 to 0.1, and it is adaptively adjusted according to the size of the patch. For the severely curved sections, a cubic polynomial was used for fitting: , ( , , (Solved using boundary smoothing constraints); For the smooth sections, keep the original boundary pixels unchanged; S4-3 stitches the fitted boundary pixels with the original smooth segment pixels to obtain the corrected patch boundary.

10. A method for rapid selection of remote sensing images for routine monitoring of natural resources according to claim 9, characterized in that, The accuracy verification in S4 uses a hierarchical accuracy weighted calculation method, including: S4-4 divides the changed patch dataset into three sub-datasets according to the semantic annotation in S1: natural change sub-dataset Repairing changed subsets of data Abnormal change subset Calculate the confusion matrix for each subset of data to obtain the overall accuracy within each class. , , The value range is 0-1; S4-5 sets the weight coefficients for the three types of subsets based on the weight requirements of the monitoring business. , , Abnormal changes have the highest weight and require close monitoring; natural changes have the lowest weight because misjudgment has a smaller impact, and are weighted using a weighted formula. Calculate the overall accuracy; S4-6 uses the same weighting coefficients to calculate the stratified Kappa coefficients. , These are the Kappa coefficients for the three subsets of data; S4-7if and If the accuracy meets the standard, then the process is considered to be satisfactory; otherwise, return to S3 for reprocessing.

Images