Method for marking quality of pixels of domestic satellite multispectral image

Abstract

This application relates to the field of remote sensing image processing technology, and provides a pixel-by-pixel quality labeling method for multispectral images from domestic satellites. Addressing the issues of unstable quality labeling accuracy and difficulty in segmenting clouds and snow in multispectral images from domestic satellites, the algorithm in this application, based on thresholding combined with guided filtering and connected region pixel matching, divides the land surface into three categories—water areas, flat areas, and mountainous areas—using numerical elevation data for quality labeling. It corrects the cloud detection threshold using heterogeneous low-resolution surface reflectance data products, corrects false detections of snow-capped mountains in mountainous areas using heterogeneous low-resolution snow cover data products and digital elevation data, and corrects false detections of snowfall in flat areas using heterogeneous low-resolution snow cover data products. Finally, it obtains a quality-labeled image containing categories for clouds, snow, clear land surface, water bodies, and fill values. The key steps involved in this application are implemented using mature algorithms, exhibiting high stability.

Description

Technical Field

[0001] This application relates to the field of remote sensing image processing technology, specifically to a pixel-by-pixel quality marking method for multispectral images from domestic satellites. Background Technology

[0002] The massive amounts of multispectral image data acquired by remote sensing satellites are processed through complex, customized algorithms to produce standardized data products that are convenient for users in various fields. In recent years, the quality of standard remote sensing satellite image data products has greatly improved, with enhanced quantitative accuracy. Pixel-by-pixel quality-labeled data products have emerged, enabling individual labeling of types such as clouds, cloud shadows, land, water bodies, and snow cover. This allows data users to filter out interfering pixels according to their application scenarios, such as filtering clouds and cloud shadows when studying surface changes and classification. Currently, commonly used data sources internationally, including MODIS, Landsat TM / ETM+ / OLI, and Sentinel, all include quality-labeled data products. These quality-labeled data products are typically provided to data users as components of standard data products, either as independent bands or images.

[0003] The core of the quality labeling algorithm in standard data products of remote sensing satellite multispectral images is cloud and cloud shadow detection. Data from the International Satellite Cloud Climatology Project (ISCCP) shows that over 50% of global satellite remote sensing data is covered by clouds. Cloud / cloud shadow detection and labeling are crucial prerequisites for efficient remote sensing applications. Furthermore, labeling categories such as water bodies and snow cover is also highly valuable for subsequent remote sensing applications. Thresholding is currently the most prevalent cloud detection method, utilizing the differences in reflectance between clouds and typical ground features in the visible and near-infrared bands, and temperature differences in the thermal infrared band, to identify clouds. It has achieved relatively high accuracy on foreign satellite data. Typical fixed-threshold cloud detection algorithms include the ISCCP method, the CLAVR method, and the APOLLO method. Represented by the US Landsat series of satellites, the QA quality-labeled band data included in its standard data products Collection 1 / 2 uses CFMAS. [1] The cloud detection algorithm comprehensively utilizes visible light, near-infrared, short-wave infrared, and thermal infrared bands, along with threshold values, to distinguish clouds, shadows beneath clouds, bright ground surfaces, snow, etc., resulting in improved accuracy compared to the original ACCA algorithm.

[0004] China is currently a major player in remote sensing satellite launches, with medium- and high-resolution multispectral sensors serving as the primary imaging payload. In recent years, many satellite data standard products have begun to include pixel-by-pixel quality labeling data. Developing quality labeling algorithms for the latest domestic satellite data has become a task and challenge for many Chinese scholars. Compared to foreign satellite data, the main gap in domestic satellite multispectral data lies in the limited number of bands and the low level of quantification. Domestic satellite multispectral images are mostly four-band data ranging from visible light to near-infrared. Due to the lack of a mid-infrared band, quality labeling algorithms cannot utilize the low-temperature characteristics of high clouds in the infrared band. This makes it difficult to effectively distinguish between clouds and the spectral similarities of bright surfaces such as snow and deserts, becoming a major source of error in domestic satellite quality labeling algorithms.

[0005] Currently, there are three main approaches to improving the accuracy of domestically developed multispectral satellite quality labeling algorithms: one is to combine image processing algorithms to improve accuracy, such as the approach taken by Shen Huanfeng et al. at Wuhan University. [2] By using guided filtering to correct cloud edges, the accuracy of edge regions in WFV multispectral images from the domestic GF-1 / 2 satellites was improved. Secondly, accuracy was further enhanced by incorporating multi-source auxiliary data, such as Lin et al. [3] A dynamic threshold cloud detection algorithm (UDTCDA) supported by a priori surface reflectance database is proposed. Thirdly, deep learning algorithms are used to train the model through a large number of data samples, enabling the model to directly distinguish targets such as clouds and snow in RGB color images, achieving a level of performance similar to human interpretation. For example, the Cloud-AttU algorithm based on the U-Net network... [4] Cloud detection methods rely on unique features for identification. These studies provide valuable insights for domestically developed multispectral satellite quality labeling standards.

[0006] References:

[0007] [1] Zhu, Zhe, and Curtis E. Woodcock. 2012. "Object-based cloud and cloudshadow detection in Landsat imagery." Remote Sensing of Environment 118:83-94.doi:10.1016 / j.rse.2011.10.028.

[0008] [2] Zhiwei Li, Huanfeng Shen, Huifang Li, Guisong Xia, Paolo Gamba, Liangpei Zhang. Multi-feature combined cloud and cloud shadow detection in GaoFen-1 wide field of view imagery[J]. Remote Sensing of Environment 191(2017)342-358.

[0009] [3]Lin sun, jing wei, jian wang, et al.A Universal Dynamic ThresholdCloud Detection Algorithm(UDTCDA) supported by a prior surface reflectancedatabase[J].Journal of Geophysical Research, D.Atmospheres:JGR, 2016, 121(12):7172-7196.DOI:10.1002 / 2015JD024722.

[0010] [4] Guo, Yanan, Xiaoqun Cao, Bainian Liu, and Mei Gao. 2020. "CloudDetection for Satellite Imagery Using Attention-Based U-Net ConvolutionalNeural Network." Symmetry 12(6):1056.doi:10.3390 / sym12061056. Summary of the Invention

[0011] This application provides a pixel-by-pixel quality marking method for domestic satellite multispectral images. The purpose is to provide an automatic pixel-by-pixel quality marking algorithm for data preprocessing applications of domestic satellite multispectral images, especially for the production of multispectral image standard data products.

[0012] This application provides a pixel-by-pixel quality marking method for domestically produced satellite multispectral images, including:

[0013] Reference data preparation: Based on the input multispectral image, prepare low-resolution satellite standard surface reflectance products, snow cover products, and DEM elevation images that cover the geographical area and are close to the imaging date, and perform cropping and reprojection to obtain reference data with the same projection, resolution, and pixel size as the input image.

[0014] The surface was divided into three categories—water area, flat area, and mountainous area—using DEM elevation images. Different quality labels were applied to the three categories to generate single-band quality label mask images.

[0015] For cloud detection, different fixed thresholds are used to detect thick clouds with high confidence in three categories: water areas, flat areas, and mountainous areas. The clouds are then marked in the quality labeling mask image.

[0016] Dynamic threshold cloud correction: For flat and mountainous areas, a linear function is constructed between the apparent reflectance of the input image and each band of the low-resolution surface reflectance data in the reference data. If the apparent reflectance of the cloud pixel in the quality marker mask is less than the maximum value of the sharp surface apparent reflectance calculated by the linear function in each band, then the cloud pixel is corrected to a sharp surface. Large-scale guided filtering is used to correct the cloud edge for the remaining cloud pixels.

[0017] Step 500, Snow Mountain Detection: For mountainous areas, use the low-resolution snow cover data in the reference data and the DEM elevation image to determine whether there are snow mountains. If there are, perform snow mountain simulation to obtain a simulated snow mountain mask, and use the simulated snow mountain mask to correct the false cloud areas in the quality marker mask of the mountainous area to snow.

[0018] Shadow detection and cloud shadow matching: For flat and mountainous areas, a fixed threshold is used to detect shadow pixels with high confidence, and a guided filter is used to correct the shadow edges. The azimuth and distance range of the shadow under the cloud relative to the cloud are calculated using imaging geometry. The correspondence between the cloud and the shadow under the cloud is determined by matching the maximum number of pixels in the connected area between the cloud and the shadow under the cloud, and pixels in the connected area of ​​the shadow without matching relationship are filtered out.

[0019] For snowfall detection, in flat areas, and for mid-to-high latitude winter snowfall areas, the low-resolution snow cover data in the reference data is used to correct the connected areas of cloud pixels that do not match shadows in the quality marker mask to snow.

[0020] Mask integration yields a quality-labeled mask image that includes clouds, snow, clear ground, water bodies, and fill value categories.

[0021] In one embodiment, the method of dividing the land surface into three categories—water areas, flat areas, and mountainous areas—using DEM elevation images includes:

[0022] Flat areas and mountainous areas are distinguished by combining the flood filling algorithm with the statistics of slope and height; different quality labels are used for the three categories. In the byte-type quality label mask image with a value range of 0 to 255, the three categories are labeled in the tens place, and the other categories are labeled in the units place.

[0023] In one embodiment, for flat and mountainous areas, constructing a linear function between the apparent reflectance of the input image and the low-resolution surface reflectance data in the reference data for each band involves fitting the linear functional relationship between apparent reflectance and surface reflectance using the least squares method. This results in a dynamic threshold cloud detection model for four channels: visible light and near-infrared. The data point set used for fitting satisfies the condition that a low-resolution surface reflectance pixel corresponds to a pixel region of a clear surface in the high-resolution input image (flat or mountainous area), and that this pixel region is an approximately homogeneous surface. If the apparent reflectance of a cloud region pixel in the quality marker mask is less than the maximum value of the apparent reflectance of a clear surface calculated by the linear function in each band, then that pixel is corrected to be a clear surface.

[0024] In one embodiment, the step of correcting the connected regions of cloud pixels in the quality marker mask that do not match shadows to snow using low-resolution snow cover data in the reference data involves first determining whether there are cases where low-resolution snow cover data snow pixels correspond to cloud pixels in the flat area of ​​the quality marker mask. If such cases exist, and the connected regions where the cloud pixels are located do not match shadows, then the connected regions of the cloud pixels are corrected to snow.

[0025] The pixel-by-pixel quality labeling method for domestic satellite multispectral images provided in this application addresses the problems of unstable quality labeling accuracy and difficulty in separating clouds and snow in domestic satellite multispectral images. Based on thresholding combined with guided filtering and connected region pixel matching algorithms, the algorithm uses numerical elevation to divide the land surface into three categories: water areas, flat areas, and mountainous areas, and performs quality labeling processing for each category. It uses heterogeneous low-resolution surface reflectance data products to correct cloud detection thresholds, heterogeneous low-resolution snow cover data products and digital elevation to correct false detections of snow-capped mountains in mountainous areas, and heterogeneous low-resolution snow cover data products to correct false detections of snowfall in flat areas. Finally, it obtains a quality-labeled image containing categories for clouds, snow, clear land surface, water bodies, and fill values. The key steps involved in this invention are implemented using mature algorithms and have high stability.

[0026] Compared with existing technologies, this application has the following characteristics: This invention provides a solution for a quality labeling algorithm for multispectral images from domestic satellites. The algorithm is fully automated, requiring no human-computer interaction throughout the process; users only need to perform a simple check on the final detection results. Key steps utilize mature algorithms, exhibiting high stability and applicability. It provides crucial technical support for the automated production of high-precision quality-labeled data products from massive amounts of domestic satellite data. This technology is a necessary data preprocessing step in the research and production of standard products for domestic satellite data. By introducing multi-source reference data, it addresses the instability of fixed-threshold cloud detection accuracy and the spectral inseparability of clouds and snow in the visible and near-infrared bands. It provides a feasible technical solution for producing quality-labeled data products from domestic satellite data. Attached Figure Description

[0027] To more clearly illustrate the technical solutions in this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0028] Figure 1 This is a flowchart of the domestic satellite multispectral image quality marking algorithm provided in the embodiments of this application;

[0029] Figure 2 This is a schematic diagram of the large-scale guided filtering edge correction method for pixel-by-pixel quality marking of domestic satellite multispectral images provided in the embodiments of this application;

[0030] Figure 3 This is a schematic diagram of pixel matching between clouds and shadow areas below clouds in the pixel-by-pixel quality marking method for domestic satellite multispectral images provided in this application embodiment;

[0031] Figure 4 This is an example diagram of the algorithm for processing Gaofen-1 satellite data using the pixel-by-pixel quality marking method for multispectral images of domestic satellites provided in this application embodiment. Detailed Implementation

[0032] To make the objectives, technical solutions, and advantages of this application clearer, the technical solutions of this application will be clearly and completely described below with reference to the accompanying drawings of the embodiments. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0033] This technology is based on two limitations of current domestic satellite multispectral remote sensing image data quality labeling algorithms: 1) Fixed threshold methods have limited cloud detection accuracy and are even less effective at defining precise cloud boundaries. 2) The effective bands for cloud detection are lacking, only including visible and near-infrared bands, failing to utilize the low-temperature characteristics of thick clouds in the infrared band, leading to numerous false detections of snow and bright surfaces. To address these two issues, and drawing on existing research, this technology designs an automated domestic satellite multispectral image quality labeling algorithm. It dynamically corrects the cloud detection threshold using heterogeneous reference data and distinguishes between false detections of snow-capped mountains and snowfall, thus meeting the need for automatically generating high-precision pixel-by-pixel quality labeling products from massive amounts of domestic satellite data.

[0034] The pixel-by-pixel quality marking method for domestic satellite multispectral images provided by the present invention will be described in detail below with reference to embodiments.

[0035] Figure 1 A flowchart of the domestic satellite multispectral image quality labeling algorithm provided in this application embodiment. (Refer to...) Figure 1 This application provides a pixel-by-pixel quality marking method for domestically produced satellite multispectral images, which may include:

[0036] Reference data preparation: Based on the input multispectral image, prepare low-resolution satellite standard surface reflectance products, snow cover products, and DEM (Digital Elevation Model) elevation images that cover the geographical area and are close to the imaging date. Then, crop and reproject to obtain reference data with the same projection, resolution, and pixel size as the input image.

[0037] The surface was divided into three categories—water area, flat area, and mountainous area—using DEM elevation images. Different quality labels were applied to the three categories to generate single-band quality label mask images.

[0038] For cloud detection, different fixed thresholds are used to detect thick clouds with high confidence in three categories: water areas, flat areas, and mountainous areas. The clouds are then marked in the quality labeling mask image.

[0039] Dynamic threshold cloud correction: For flat and mountainous areas, a linear function is constructed between the apparent reflectance of the input image and each band of the low-resolution surface reflectance data in the reference data. If the apparent reflectance of the cloud pixel in the quality marker mask is less than the maximum value of the sharp surface apparent reflectance calculated by the linear function in each band, then the cloud pixel is corrected to a sharp surface. Large-scale guided filtering is used to correct the cloud edge for the remaining cloud pixels.

[0040] For snow mountain detection, in mountainous areas, the low-resolution snow cover data in the reference data and the DEM elevation image are used to determine whether snow mountains exist. If they exist, a snow mountain simulation is performed to obtain a simulated snow mountain mask. The simulated snow mountain mask is then used to correct the false cloud areas in the quality marker mask of the mountainous area to snow.

[0041] Shadow detection and cloud shadow matching: For flat and mountainous areas, a fixed threshold is used to detect shadow pixels with high confidence, and a guided filter is used to correct the shadow edges. The azimuth and distance range of the shadow under the cloud relative to the cloud are calculated using imaging geometry. The correspondence between the cloud and the shadow under the cloud is determined by matching the maximum number of pixels in the connected area between the cloud and the shadow under the cloud, and pixels in the connected area of ​​the shadow without matching relationship are filtered out.

[0042] For snowfall detection, in flat areas, and for mid-to-high latitude winter snowfall areas, the low-resolution snow cover data in the reference data is used to correct the connected areas of cloud pixels that do not match shadows in the quality marker mask to snow.

[0043] Mask integration yields a quality-labeled mask image that includes clouds, snow, clear ground, water bodies, and fill value categories.

[0044] In one embodiment, the method of dividing the land surface into three categories—water areas, flat areas, and mountainous areas—using DEM elevation images includes:

[0045] Flat areas and mountainous areas are distinguished by combining the flood filling algorithm with the statistics of slope and height; different quality labels are used for the three categories. In the byte-type quality label mask image with a value range of 0 to 255, the three categories are labeled in the tens place, and the other categories are labeled in the units place.

[0046] In one embodiment, for flat and mountainous areas, constructing a linear function between the apparent reflectance of the input image and the low-resolution surface reflectance data in the reference data for each band involves fitting the linear functional relationship between apparent reflectance and surface reflectance using the least squares method. This results in a dynamic threshold cloud detection model for four channels: visible light and near-infrared. The data point set used for fitting satisfies the condition that a low-resolution surface reflectance pixel corresponds to a pixel region of a clear surface in the high-resolution input image (flat or mountainous area), and that this pixel region is an approximately homogeneous surface. If the apparent reflectance of a cloud region pixel in the quality marker mask is less than the maximum value of the apparent reflectance of a clear surface calculated by the linear function in each band, then that pixel is corrected to be a clear surface.

[0047] In one embodiment, the step of correcting the connected regions of cloud pixels in the quality marker mask that do not match shadows to snow using low-resolution snow cover data in the reference data involves first determining whether there are cases where low-resolution snow cover data snow pixels correspond to cloud pixels in the flat area of ​​the quality marker mask. If such cases exist, and the connected regions where the cloud pixels are located do not match shadows, then the connected regions of the cloud pixels are corrected to snow.

[0048] This application includes the following:

[0049] Input data: Domestic multispectral images, specifically image data acquired by medium-to-high resolution multispectral sensors onboard domestic military / civilian satellites, such as the GF series, HJ series, ZY series, and JB series satellites. The acquired multispectral images have a spatial resolution of no less than 30 meters and contain only four bands from visible light to near-infrared. The input data has been converted from raw DN (Digital Number, remote sensing image pixel brightness value) to apparent reflectance using calibration parameters.

[0050] Processing steps: The algorithm automatically queries and downloads reference data. Based on the geographical range and imaging time of the input multispectral image, it automatically downloads low-resolution satellite standard surface reflectance and snow cover products with similar imaging dates and geographical coverage online. It extracts the corresponding surface reflectance band and snow cover band, and then crops and reprojects them to obtain a Byte-type image with the same projection, resolution, and pixel size as the input image.

[0051] Low-resolution satellite standard surface reflectance products are currently available using the 500-meter resolution MOD09A1 product from the MODIS series satellites. This data source can be automatically queried and downloaded online after account registration, utilizing the geographic range and imaging time code of the input image. This product is based on high observation coverage, low viewing angle, no clouds or cloud shadows, and aerosol concentration for precise positioning and geometric correction. It selects the best Level 2 Grid for each product pixel within 8 days for observation, resulting in high overall accuracy. The MOD09A1 product is used in the blue light band (0.459–0.479 μm), green light band (0.545–0.565 μm), red light band (0.62–0.67 μm), and near-infrared band (0.841–0.876 μm).

[0052] Low-resolution satellite standard snow cover products can currently use the 500-meter resolution standard snow cover products from the MODIS series satellites. This is achieved through automatic online querying and downloading, downloading either 8-day composite MOD10A2 data or daily MOD10A1 data, using the Fractional Snow Cover (FSC) band. The MOD10A1 snow cover product is generated based on multi-band sensor imagery from MODIS series satellites. Automatic snow cover mapping is performed using the Normalized Snow Index (NSDI) calculated from the reflectance coefficients of the satellite in MODIS bands 4 (0.545–0.565 μm) and 6 (1.628–1.652 μm), combined with threshold testing and MODIS cloud mask data, generating a binary daily-scale snow cover product (“Snow Cover Daily Tile”). MOD10A2 is an 8-day synthetic data product obtained through further integration and processing. It utilizes the motion characteristics of clouds to filter out some cloud coverage, resulting in higher accuracy and quality. However, there are cases where some product data cannot be generated or retrieved.

[0053] The recommended spatial resolution for DEM digital elevation is no less than 90 meters. Currently, global SRTM data with a spatial resolution of 30 meters can be used. Unlike MODIS data products, SRTM data is not updated daily and requires online querying and downloading when needed. All SRTM data can be pre-downloaded to the local hard drive for direct querying and retrieval. Based on the geographical extent of the input multispectral image, the DEM digital elevation is cropped and reprojected to obtain a single-band DEM elevation image with the same projection, resolution, and pixel size as the input image.

[0054] Processing steps: Elevation classification. Based on the DEM elevation image, the land surface is divided into three categories: water areas, flat areas, and mountainous areas. The purpose is to use different quality labeling algorithms for different land surface categories to improve the overall quality labeling accuracy. It should be noted that the above three categories are a simple and coarse division, which is different from precise remote sensing classification. The division of the three categories only serves the quality labeling algorithm process and is not reflected in the final quality-labeled data product.

[0055] Water body classification is determined based on the data products of the DEM elevation image. For example, the -32768 fill value in the SRTM data is marked as a water body. This water body category only includes large areas of deep sea and lakes, generally areas where SAR cannot measure topography, and does not include small lakes and rivers. This is simple and applicable for threshold-based cloud and shadow detection. It should be noted that the final quality labeling data product includes water body labels using other reference data, combined with specialized fusion and water body detection algorithms. The reference data for water body labeling uses Global Surface Water Occurrence (GSWO) data, which can also be downloaded automatically from the network.

[0056] Flat areas and mountainous areas are distinguished using the Flood Fill algorithm. Flood Fill is a classic image segmentation algorithm, named for its analogy to a flood spreading from one region to all reachable areas. In GNU Go and Minesweeper, Flood Fill is used to calculate the areas that need to be cleared. Flood Fill can be constructed in various ways, and many implementation code resources exist. Most algorithms explicitly or implicitly use queue or stack data structures. The Flood Fill algorithm has two main parameters: the starting point and the target height. The starting point and target height are determined by statistically analyzing the slope and height of a DEM (Digital Elevation Model) image. The pixel values ​​in the DEM image represent the height, and the slope is calculated from the height difference of a sliding window and the pixel size of the sliding window. A fixed slope threshold is given as the boundary between flat and steep areas. A histogram of all pixel values ​​below this threshold determines the maximum starting point, and a histogram of all pixel values ​​above this threshold determines the maximum target height.

[0057] The three categories use different quality labels in the quality label mask image. In the Byte-type quality label mask image with a value range of 0 to 255, the three categories are labeled in the tens place, for example, water is labeled 10, flat areas are labeled 20, and mountainous areas are labeled 30. Other categories are labeled in the units place, for example, fill value is labeled 0, clear surface is labeled 1, snow is labeled 4, and cloud is labeled 5. For example, the value 15 represents clouds in water areas.

[0058] The core of this algorithm's accuracy improvement lies in classifying regions into three categories: water areas, flat areas, and mountainous areas, and then separately labeling them for quality. Using only DEM elevation data for these three categories is a simple and effective engineering technique. Previous cloud detection algorithms, such as the SFM algorithm designed for the Gaofen-1 satellite, used simple spectral thresholds to distinguish between land and water. However, due to the limitations of multispectral image band settings on domestic satellites, accuracy became unstable. Subsequent algorithms utilized multi-source data and designed complex processing steps to construct reference datasets. While these methods achieved high quantitative accuracy, they were difficult to apply to automated global data processing workflows. Furthermore, experiments show that for quality labeling algorithms centered on cloud detection, the accuracy requirement for category classification is not necessarily higher with more precise classification, but rather depends more on how the overall brightness, texture, and topography of the underlying land surface differentiates cloud thresholds and shadow features. This algorithm, by using DEM elevation data for classification, typically achieves better accuracy.

[0059] Processing steps: Cloud detection. For the three categories of water areas, flat areas, and mountainous areas, different fixed thresholds are used to detect thick clouds with high confidence. For water areas, large-scale guided filtering is used to correct the cloud area edges, and the clouds are marked in the quality labeling mask.

[0060] The input multispectral image is converted from raw DN values ​​to apparent reflectance. The fixed threshold method for cloud detection sets a threshold for each band; when the pixel value at a given location exceeds the threshold for all bands, that pixel is marked as a cloud. Different thresholds are used for three categories—water areas, flat areas, and mountainous areas—to detect thick clouds with high confidence. The specific thresholds are determined based on satellite data and experimental accumulation and are set as algorithm configuration items that can be manually modified at any time. To obtain reasonable accuracy for data from different times and geographical areas, the fixed threshold is usually set strictly, ensuring that the detection results primarily show thick clouds.

[0061] The cloud detection uses both the HOT index and the VBR index, with the detection results primarily based on thick clouds.

[0062] The formula for calculating the HOT index is:

[0063] HOT = B1 - 0.5 × B3

[0064] The formula for calculating the VBR index is:

[0065]

[0066] B1, B2, and B3 represent the blue, green, and red band data, respectively. Different thresholds are used for the cloud detection index for the three categories: water areas, flat areas, and mountainous areas. For example, the threshold for flat areas is HOT>0.2 && VBR>0.7, while the thresholds for mountainous and water areas decrease progressively. Clouds are marked in the quality marker mask, with water cloud areas marked as 15 and surface cloud areas marked as 25 or 35.

[0067] Processing Steps: Dynamic Threshold Cloud Correction. Flat and mountainous areas are treated as a single region. Clouds detected by a fixed threshold in the quality-marked mask within this region have low accuracy and may contain numerous false detections. To refine the cloud detection threshold for the current data, low-resolution surface reflectance data from similar imaging dates are used as reference data to find a more accurate cloud detection threshold. Specifically, the steps are as follows: First, a linear function is constructed between the apparent reflectance of the input image and the clear surface in the low-resolution surface reflectance data from the reference data for each band. Then, based on whether the apparent reflectance of cloud pixels in the quality-marked mask is greater than the linear function in each band, the maximum value of the apparent reflectance of the clear surface is calculated to determine if the pixel is a corrected cloud. Finally, large-scale guided filtering is used to correct the cloud area edges.

[0068] Construct a linear function yi = ax for each band between the apparent reflectance of the input image and the low-resolution surface reflectance data of the reference data (clear surface). i +b, where x i Corresponding to low-resolution surface reflectance, y i The apparent reflectance corresponds to the input image, i = 1, 2, 3… corresponding to band numbers, and a and b are the linear variation coefficients to be solved. Each band is different and needs to be calculated separately. If the reference data point set {x} is determined… i} and the corresponding input image point set {y i Then, the least squares method can be used to fit the linear functional relationship between apparent reflectance and surface reflectance, and a dynamic threshold cloud detection model with four channels of visible light and near-infrared light can be constructed.

[0069] The data point set {x} used for fitting i y i A point set {x} satisfies the condition that a low-resolution surface reflectance pixel corresponds to a pixel region of clear surface in both flat and mountainous areas of a high-resolution input image, and that this pixel region is an approximately homogeneous surface. i y iThe specific method for obtaining the data is illustrated using the 16-meter resolution domestic Gaofen-1 (GF-1) and the 500-meter resolution MODIS MOD09A1 as examples. First, a possible point set is selected in MOD09A1. A quality marker mask is used to select all pixels in flat and mountainous areas. Cloud pixels and low-quality pixels are filtered out based on the quality marker bands in MOD09A1. Points with large pixel data differences (generally boundary points) are filtered out based on the difference in pixel values ​​between their eight neighbors. The remaining points constitute the possible point set in MOD09A1. Then, for all possible point sets in MOD09A1, each pixel at the 500-meter resolution is mapped to multiple pixels within a rectangular frame at the 16-meter resolution of GF-1, based on their geographic coordinates. All pixels within the rectangular frame are counted. If no cloud pixels are included and the pixel difference value is below a threshold, the point is selected into the point set {x}. i ,yi}, and take the midpoint of the rectangle's pixels as y i Value. The GF-1 multispectral image contains four bands. The selected point set uses only the fourth near-infrared band, which corresponds to a band with a similar center band length in MOD09A1.

[0070] Fitting the point set {x} using the least squares method i y i If a linear function fitting fails for a certain band, or if the point set {x}... i y i If the mean square error of the fitted linear function is greater than the threshold, dynamic threshold correction is abandoned, and the detection results with a fixed threshold are directly used for guided filtering to correct the cloud edges. If the fitting is successful, the apparent reflectance cloud threshold of GF-1 is calculated using the fitted linear function and the surface reflectance band of MOD09A1. The specific steps are: first, all pixels (including cloud pixels) in the surface reflectance band of MOD09A1 are used to calculate the apparent reflectance using the fitted linear function; then, the threshold between the surface and cloud in the obtained apparent reflectance is statistically analyzed, and this threshold is used as the dynamic threshold for cloud detection correction. The calculation of the dynamic threshold for cloud detection uses all four bands contained in the GF-1 multispectral image, corresponding to four bands with similar center band lengths in MOD09A1.

[0071] The quality marker mask is corrected using a dynamic threshold based on cloud detection. If the apparent reflectance of a pixel in the cloud area of ​​the quality marker mask is less than the maximum value of the apparent reflectance of a sharp surface calculated by a linear function in all bands, then that pixel is corrected to a sharp surface. The corrected result uses a large-scale guided filter to correct the cloud area edges. Specifically, the quality marker mask image is used as input, the blue light band is used as the guide map, and the guided filter algorithm adds thin clouds that were missed at the cloud area edges to the cloud mask. Large-scale refers to a filter radius pixel size greater than 500 pixels, which is intended to accommodate the correction of missed detections in large-area cloud areas. In the guided filter algorithm, the average band of the three visible light bands is used as the guide map L(x,y), the cloud mask is used as the filter input image V(x,y), and the output image is denoted as... The filtering result at pixel (x, y) is expressed as a weighted average:

[0072]

[0073] Assume the guided filter has a guided image L(x,y) and a filtered output. The relationship is a locally linear model:

[0074]

[0075] By minimizing the following window ω k Cost function:

[0076]

[0077] The local linear coefficients a are obtained k ,b k The value of .

[0078] The main computational burden of the guided filter is the median filter. To accommodate the transition region between large areas of thick clouds and the ground surface, i.e., cloud boundaries or thin cloud areas, a large filtering radius of 500 pixels is used. Increasing the filtering radius leads to a sharp increase in the computational efficiency of a typical median filter; therefore, a fast filtering algorithm independent of the radius is adopted using the Boxfilter. This reduces the complexity of operations such as summation, mean, and variance (O(MN)) to approximately O(1). The output image of the guided filter... For corrected cloud area masking.

[0079] Figure 2 This is a schematic diagram of large-scale guided filtering edge correction in cloud detection, showing that guided filtering of images ( Figure 2 .c), the resulting cloud mask ( Figure 2 .d) compared to before ( Figure 2 b) More consistent with the original image Figure 2 The actual cloud region boundary in .a).

[0080] Processing steps: Snow mountain detection. Due to the lack of effective bands for distinguishing between clouds and snow in domestic satellite multispectral data, threshold-based cloud detection is prone to misdetecting snow as clouds. Snow mountains, as a type of snow cover, can be distinguished through snow mountain simulation because they are related to terrain elevation. For mountainous areas, low-resolution snow cover data from the reference data and DEM elevation images are used to determine the presence of snow mountains. If present, a snow mountain simulation is performed to obtain a simulated snow mountain mask. This simulated snow mountain mask is then used to correct misdetected cloud areas in the quality-marked mask for snow cover.

[0081] The presence of snow-capped mountains is determined by using low-resolution snow cover data from the reference data and DEM elevation images. If the mountainous area in the DEM elevation image contains low-resolution snow cover pixels and the proportion of snow-capped mountain pixels exceeds 2%, then a snow-capped mountain simulation is performed.

[0082] The first step in simulating a snow mountain is to calculate the average elevation value of the mountain area pixels in the DEM image corresponding to all snow pixels in the snow mask, which is then used as the snow line value. Pixels in the DEM image with a value greater than the snow line value are marked as snow mountains in the simulated snow mountain mask. The simulated snow mountain mask is a byte-type single-band image with the same projection, resolution, and pixel size as the input image.

[0083] Using a simulated snow mountain mask to correct falsely detected cloud areas in the mountainous area of ​​the quality marker mask to snow involves performing morphological matching between all cloud pixel connected regions in the mountainous area and the simulated snow mountain pixel connected regions in the simulated snow mountain mask. If the similarity match exceeds a threshold, the cloud pixel connected regions are corrected to snow.

[0084] Processing steps: Shadow detection and cloud shadow matching. Under-cloud shadows usually do not contain surface information or only contain weak surface information. The quality labeling algorithm correctly identifies under-cloud shadows, which is meaningful for the subsequent application of labeled data products.

[0085] Shadow detection uses a thresholding approach combined with image processing algorithms to mark all possible shadow areas in flat and mountainous regions, including shadows under clouds and other dark surfaces. First, high-confidence shadow pixels are detected using a fixed threshold. Then, guided filtering is used to correct shadow edges. The fixed threshold for shadow pixel detection uses a near-infrared band threshold, NDVI index, and VBR index to set a fixed threshold. The NDVI index calculation formula is as follows:

[0086]

[0087] B3 and B4 represent the red band and near-infrared band data, respectively. For example, the threshold for coarse shadow detection is NDVI < 0.15 && VBR > 0.7 && B4 < 0.15.

[0088] When using guided filtering to correct shadow edges, unlike cloud edge correction, shadow edge correction does not result in missed detections in large radius areas. Therefore, a maximum radius of 50 pixels for guided filtering is sufficient. Additionally, floodfill is a commonly used algorithm for shadow area edge correction. However, because the algorithm expands unrestricted before the seed point reaches its highest "water level," guided filtering was chosen to correct for missed edge detections when shadow threshold detection accuracy is normal.

[0089] Cloud shadow matching utilizes morphological image processing algorithms to filter out under-cloud shadows from cloud shadow detection results. The azimuth angle of the under-cloud shadow relative to the cloud is calculated using imaging geometry. Solar zenith angle, solar azimuth angle, observed zenith angle, and observed azimuth angle can be obtained from satellite data auxiliary files. The formula for calculating the azimuth angle of the under-cloud shadow relative to the cloud area is:

[0090]

[0091] Where α, β, ω, and μ are the solar zenith angle, solar azimuth angle, observed zenith angle, and observed azimuth angle, respectively.

[0092] The distance range of the under-cloud shadow relative to the cloud is determined empirically. For example, in a 16-meter resolution multispectral image of the GF-1 satellite, the distance range of the cloud shadow is set to 0–300 pixels. For each cloud pixel connected region, the under-cloud shadow is moved 0–300 pixels in the azimuth direction relative to the cloud, moving one pixel at a time, and the number of currently matching pixels is calculated. The connected region of the shadow with the maximum number of matching pixels is marked as the under-cloud shadow. Figure 3 This is a schematic diagram of pixel matching between clouds and the shadow area beneath them. The actual algorithm was optimized to improve the computational efficiency of cloud-shadow matching, including resampling the quality marker mask to half the size of the original image, matching only connected regions with more than 50 cloud pixels, and a single matching step size of 4 pixels.

[0093] For large water areas, this algorithm does not perform under-cloud shadow detection. Considering that the overall brightness of water is relatively low compared to the ground surface, the accuracy of threshold-based shadow detection is unstable, the matching of clouds and under-cloud shadows is easily affected by water texture, and the application scenarios of water shadows in quality-labeled data are limited, therefore, water areas in the quality-labeled mask do not include under-cloud shadow markers.

[0094] Processing steps: Snowfall detection, primarily targeting data from mid-to-high latitude winter snowfall areas, determines whether large-scale snowfall exists in flat areas. This is done by using low-resolution snow cover data from a reference dataset to assess the presence of large-scale snowfall in these areas. The criterion is that, after excluding cloud pixels, the proportion of snow pixels in the image is greater than that of the clear ground surface. If large-scale snowfall exists in a flat area, connected regions of cloud pixels that do not match shadows are corrected to snow cover, while connected regions of cloud pixels that match shadows below clouds are considered clouds above the snowfall area.

[0095] Processing steps: Mask integration is to integrate the quality-labeled mask data to conform to the specifications of the quality-labeled data product. After regional integration into three categories—water area, flat area, and mountain area—the final quality-labeled mask includes cloud, snow, clear surface, water body, and fill value categories.

[0096] Output data: Quality marker mask, which is byte-type single-band compressed TIFF format image data with the same pixel size as the input image. It includes cloud, snow, clear ground, water body, and fill value categories. Different categories are marked with different values. Here, fill value is marked as 0, clear ground as 1, water body as 2, snow as 4, and cloud as 5.

[0097] It should be noted that the core of this domestic multispectral satellite quality labeling algorithm is to solve the problem of distinguishing between dynamic thresholds and clouds and snow. This algorithm is part of the production algorithm process of domestic multispectral satellite standard data products. The quality labeling mask output by this algorithm is stored on the hard disk as a component of the domestic multispectral satellite standard data products. Figure 4 This is an example image of the algorithm processing data from the Gaofen-1 satellite, specifically an example image of a 16-meter resolution WFV multispectral image. The quality marker mask produced uses a color table to distinguish the differences between different categories.

[0098] The C++ algorithm example of this invention has been implemented on the PC platform, and the effectiveness and robustness of the algorithm have been verified through experimental data in the early stage.

[0099] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application.

Claims

1. A pixel-by-pixel quality marking method for domestically produced satellite multispectral images, characterized in that, include: Reference data preparation: Based on the input multispectral image, prepare low-resolution satellite standard surface reflectance products, snow cover products, and DEM elevation images that cover the geographical area and are close to the imaging date, and perform cropping and reprojection to obtain reference data with the same projection, resolution, and pixel size as the input image. The surface was divided into three categories—water area, flat area, and mountainous area—using DEM elevation images. Different quality labels were applied to the three categories to generate single-band quality label mask images. For cloud detection, different fixed thresholds are used to detect thick clouds with high confidence in three categories: water areas, flat areas, and mountainous areas. The clouds are then marked in the quality labeling mask image. Dynamic threshold cloud correction: For flat and mountainous areas, a linear function is constructed between the apparent reflectance of the input image and each band of the low-resolution surface reflectance data in the reference data. If the apparent reflectance of the cloud pixel in the quality marker mask is less than the maximum value of the sharp surface apparent reflectance calculated by the linear function in each band, then the cloud pixel is corrected to a sharp surface. Large-scale guided filtering is used to correct the cloud edge for the remaining cloud pixels. For snow mountain detection, in mountainous areas, the low-resolution snow cover data in the reference data and the DEM elevation image are used to determine whether snow mountains exist. If they exist, a snow mountain simulation is performed to obtain a simulated snow mountain mask. The simulated snow mountain mask is then used to correct the false cloud areas in the quality marker mask of the mountainous area to snow. Shadow detection and cloud shadow matching: For flat and mountainous areas, a fixed threshold is used to detect shadow pixels with high confidence, and a guided filter is used to correct the shadow edges. The azimuth and distance range of the shadow under the cloud relative to the cloud are calculated using imaging geometry. The correspondence between the cloud and the shadow under the cloud is determined by matching the maximum number of pixels in the connected area between the cloud and the shadow under the cloud, and pixels in the connected area of ​​the shadow without matching relationship are filtered out. For snowfall detection, in flat areas, and for mid-to-high latitude winter snowfall areas, the low-resolution snow cover data in the reference data is used to correct the connected areas of cloud pixels that do not match shadows in the quality marker mask to snow. Mask integration yields a quality-labeled mask image that includes clouds, snow, clear ground, water bodies, and fill value categories.

2. The pixel-by-pixel quality marking method for domestically produced satellite multispectral images according to claim 1, characterized in that, The method of using DEM elevation images to divide the land surface into three categories: water areas, flat areas, and mountainous areas includes: Flat areas and mountainous areas are distinguished by combining the flood filling algorithm with the statistics of slope and height; different quality labels are used for the three categories. In the byte-type quality label mask image with a value range of 0 to 255, the three categories are labeled in the tens place, and the other categories are labeled in the units place.

3. The pixel-by-pixel quality marking method for domestically produced satellite multispectral images according to claim 1, characterized in that, For flat and mountainous areas, a linear function is constructed between the apparent reflectance of the input image and the low-resolution surface reflectance data in the reference data for each band. This is achieved by fitting the linear functional relationship between apparent reflectance and surface reflectance using the least squares method. A dynamic threshold cloud detection model with four channels (visible and near-infrared) is then constructed. The data point set used for fitting satisfies the condition that a low-resolution surface reflectance pixel corresponds to a pixel region of a clear surface in the high-resolution input image (flat or mountainous area), and that this pixel region is an approximately homogeneous surface. If the apparent reflectance of a cloud region pixel in the quality marker mask is less than the maximum value of the apparent reflectance of a clear surface calculated by the linear function in each band, then that pixel is corrected to be a clear surface.

4. The pixel-by-pixel quality marking method for domestically produced satellite multispectral images according to claim 2, characterized in that, The step of correcting the connected regions of cloud pixels in the quality marker mask that do not match shadows to snow using low-resolution snow cover data in the reference data involves first determining whether there are cases where low-resolution snow cover data snow pixels correspond to cloud pixels in the flat area of ​​the quality marker mask. If such cases exist, and the connected region where the cloud pixel is located does not match shadows, then the connected region of the cloud pixel is corrected to snow.

Images