A method for optimizing thin cloud thickness graph considering potential low-rank features

Abstract

The invention discloses a thin cloud thickness map optimization method considering potential low-rank features, which belongs to the technical field of remote sensing image processing. Based on remote sensing images, reflectance data of multiple bands are obtained and an initial thin cloud map is estimated. The initial thin cloud map is subjected to feature weighting processing using an optimized potential low-rank representation model of the thin cloud map to obtain a low-rank optimized cloud map. According to the reflectance data of multiple bands, a building area index map and a water body index map are calculated; and according to the building area index map and the water body index map, a non-water body building area distribution index map is calculated, a weight suppression matrix for adaptively adjusting the weight of the building area distribution index is generated, and the low-rank optimized cloud map is processed to obtain a cloud map with a reduced highlight building area, and then the cloud map is optimized to obtain a final optimized thin cloud map, so that the seams of the block edges in the thin cloud map are smoothed and the edge and gradient information of the thin cloud are enhanced. Accurate and robust technical support is provided for the estimation and correction of thin cloud thickness of remote sensing images.

Description

Technical Field

[0001] The invention relates to a thin cloud thickness map optimization method considering potential low-rank features, and belongs to the technical field of remote sensing image processing. Background Art

[0002] Optical remote sensing satellite earth observation is a complex process, which is easily affected by the atmospheric environment, thus reducing the quality of satellite remote sensing images and greatly limiting the subsequent application of remote sensing data, such as vegetation index, land use and mapping. Among the many atmospheric factors, cloud pollution is the most common and unavoidable influencing factor. Optically, clouds can be simply divided into thin clouds and thick clouds. Thick clouds refer to the inability of satellites to obtain ground information through them, resulting in the loss of ground information in remote sensing images. Thin clouds are translucent in nature, and some surface information can be observed through them from an optical point of view, resulting in a decrease in the saturation and visibility of the image, and the radiation value will change. Therefore, in order to improve the overall quality of optical remote sensing images, it is necessary to correct the optical remote sensing images contaminated by thin clouds to restore the radiation value of the remote sensing images.

[0003] In recent years, in order to correct the degraded images affected by thin clouds, many single remote sensing image thin cloud removal methods have been proposed, which can be divided into traditional mathematical model methods and deep learning-based methods. Traditional mathematical model methods can be simply divided into feature extraction-based methods, prior hypothesis-based methods, and cirrus band-assisted methods. Among them, the feature extraction-based methods have high requirements for sample selection. If the sample selection is improper, such as the surface type in the sample data is not representative, there may be problems such as incomplete thin cloud correction and spectral distortion. Among the prior hypothesis methods, the most representative ones are dark pixel subtraction (DOS) and dark channel method (DCP). Among them, the DOS method is easily affected by the image of the bright surface, making the estimated cloud map inaccurate. The DCP method is widely used in defogging in natural images and is only suitable for homogeneous thin clouds. The cirrus band-assisted method uses the cirrus band carried by the satellite, and relies on statistical relationships to solve the linear or nonlinear relationship between the reflectivity of the cirrus band and the visible-near infrared band, so as to obtain the cirrus reflectivity of each other band, thereby achieving the purpose of cirrus removal. However, in some dry scenes or low-altitude thin clouds are difficult to be detected by the cirrus band, and this method will fail. Nowadays, with the rise of deep learning, more and more deep learning-based methods have been proposed, such as those based on convolutional neural networks, generative adversarial networks, and conditional variational autoencoder networks. Although the learning-based method can obtain good thin cloud removal results, its transferability and universality are poor. When the scene switches, it is necessary to retrain the current scene data to obtain a new model. In addition, since it is difficult to obtain cloud-covered and cloud-free images at the same location and time in real situations, it is very difficult to construct its network training data set. Most current methods use simulated data sets or match data with adjacent phase images. However, the removal of thin clouds in remote sensing images involves changes in radiation values, and there will be obvious differences in radiation values ​​at different times or scenes, which greatly limits the promotion and application of deep learning models.

[0004] Whether it is the traditional mathematical model method or the method based on deep learning, they essentially use the model to accurately estimate the thickness of thin clouds and their reflectivity intensity from the images contaminated by thin clouds, and then use the image imaging model and the relationship between bands to eliminate the thin clouds in each band. Therefore, in the thin cloud elimination technology, whether the thickness of thin clouds can be accurately estimated is a crucial task, which directly affects the accuracy and correction results of remote sensing image correction. Summary of the invention

[0005] The purpose of the present invention is to provide a thin cloud thickness map optimization method taking into account potential low-rank characteristics, combining the distribution of thin clouds and building areas in remote sensing images with low-rank sparsity theory, constructing a potential low-rank representation model for thin cloud maps, and proposing a new type of building area index, which can effectively suppress the influence of building areas on thin cloud maps. The thin cloud maps can be optimized for images affected by thin clouds under different time, illumination, space and other factors, solving the problems of poor robustness and low precision of thin cloud thickness estimation results in the prior art, and providing accurate and robust technical support for thin cloud correction and thin cloud thickness estimation of remote sensing images.

[0006] In order to solve the above technical problems, the present invention is implemented by adopting the following technical solutions.

[0007] A thin cloud thickness map optimization method considering potential low-rank features, including:

[0008] Step 1: Read multispectral remote sensing thin cloud images;

[0009] Step 2: Estimate the initial thin cloud map based on the reflectance data of multiple bands in the image; the multiple bands include: coastal aerosol band, blue band, green band, red band, near infrared band and two short-wave infrared bands;

[0010] Step 3: Using the optimized thin cloud map potential low-rank representation model to perform feature weighting processing on the initial thin cloud map to obtain a low-rank optimized cloud map;

[0011] Step 4: Based on the reflectance data of multiple bands, calculate the building area index map and water body index map;

[0012] Step 5: Calculate the distribution index map of non-water body building areas based on the building area index map and the water body index map;

[0013] Step 6: Based on the non-water body built-up area distribution index map, a weight suppression matrix for adaptively adjusting the built-up area distribution index weight is generated;

[0014] Step 7: Perform weight suppression processing on the low-rank optimized cloud image according to the weight suppression matrix to obtain a cloud image with reduced highlighted building areas;

[0015] Step 8: Optimize the cloud image that reduces the highlighted building area to obtain the final optimized thin cloud image.

[0016] Optionally, in step 2, an initial thin cloud image is estimated based on reflectance data of multiple bands in the image, including:

[0017] Obtain the reflectance data of the coastal aerosol band, blue band, green band, red band, near infrared band and two short-wave infrared bands, and extract the pixels of the coastal aerosol band and blue band in the cloud-free area of ​​the image as the sample data of the linear fitting, and use other bands as the observed quantities of the linear fitting to obtain the linear fitting relationship of the sample data in the cloud-free area;

[0018] Based on the linear fitting relationship of the sample data, the M-estimation method is used to obtain the linear relationship between the reflectance of the coastal aerosol band and the blue band;

[0019] Band synthesis is performed based on the linear relationship and reflectivity between the coastal aerosol band and the blue band to obtain a synthetic band that is easily affected by thin clouds;

[0020] Processing the synthetic band susceptible to thin cloud influence based on a dark pixel search method, identifying dark pixels in a local window of the synthetic band image;

[0021] The dark pixels in the local window are resampled according to the bicubic interpolation method to generate an initial thin cloud image;

[0022] The expression for obtaining the synthetic band susceptible to thin clouds is:

[0023] (1);

[0024] In the formula, For the synthetic band, and are the reflectances of the coastal aerosol band and the blue band, and is the linear relationship coefficient obtained by using the M-estimation method for the sample data in the cloud-free area.

[0025] Optionally, the optimization method of the thin cloud map latent low-rank representation model includes:

[0026] Perform low-rank characteristic analysis on the initial thin cloud image to obtain the basic feature items and significant feature items of the initial thin cloud image; and analyze the sparse surface noise information in the initial thin cloud image;

[0027] Based on the basic feature items and significant feature items of the initial thin cloud image and the sparse surface noise information, a relationship model between low-rank thin cloud image information and sparse surface information is constructed;

[0028] According to the relationship model between low-rank thin cloud image information and sparse surface information, a convex optimization problem of features is set;

[0029] The augmented Lagrangian algorithm is used to construct the minimization optimization function for solving convex optimization problems;

[0030] Iterate multiple variables of the minimization optimization function of the convex optimization problem respectively until the convergence condition is met and stop the iteration to obtain the optimized thin cloud map potential low-rank representation model;

[0031] Among them, the expression for constructing the relationship model between low-rank thin cloud image information and sparse surface information is:

[0032] (2);

[0033] In the formula, is the initial thin cloud image, is the basic characteristic coefficient, is the basic feature item, is the significant characteristic coefficient, is a notable feature item, is the noise term, is the nuclear norm, which is used to represent The low rank property of for norm, which is used to represent Sparsity, is the regularization parameter, with a value range of [0,1];

[0034] The expression of the minimization optimization function for solving convex optimization problems is constructed as follows:

[0035] (3);

[0036] In the formula, and are the low-rank decomposition results of basic feature items and significant feature items, respectively. , and are the Lagrangian operators of noise term, basic characteristic term and significant special term, respectively. is the trace of the matrix, is the F-norm of the matrix, is the penalty parameter, .

[0037] Optionally, the method of iterating multiple variables of the minimization optimization function of the convex optimization problem respectively until the iteration is stopped when the convergence condition is satisfied includes:

[0038] Based on minimizing the variables of the optimization function, the optimization is solved, and the update sub-problems of each variable are established respectively. In the process of solving the update sub-problems of each variable, other variables are kept unchanged, and the update of the variable is stopped when the associated variable meets the convergence condition;

[0039] Among them, the variables include: low-rank decomposition results of basic feature items , the low-rank decomposition results of the significant feature items , basic characteristic coefficient , significant characteristic coefficient , noise term , the Lagrangian operator of the noise term , the Lagrangian operator of the basic characteristic term and the Lagrangian operator with significant special terms ;

[0040] The sub-problem solving process for each variable includes:

[0041] Keeping other variables unchanged, the variable The expression of the update subproblem is:

[0042] (4);

[0043] In the formula, For variables The formula for minimizing the operator, is the Lagrangian operator of the basic characteristic term, and They are in transposition relationship with each other;

[0044] variable The expressions of the update process are:

[0045]

[0046] (5);

[0047] In the formula, for Temporary variables, is the singular value decomposition method, are the three sub-matrices obtained after SVD decomposition, They are the left singular matrix, singular value matrix and right singular matrix respectively;

[0048] Using the singular value threshold method, the singular values ​​after SVD decomposition are Screen and get the number of singular values ​​that meet the conditions , in this way Update and get the expression of the number of singular values ​​that meet the conditions:

[0049]

[0050] (6);

[0051] In the formula, is the number of singular values ​​that meet the conditions, is a singular value vector;

[0052] Get variables The final update subproblem of is expressed as:

[0053]

[0054] In the formula, is a diagonal matrix;

[0055] Keeping other variables unchanged, the variable The expression of the update subproblem is:

[0056] (7);

[0057] In the formula, is the Lagrangian operator of the significant characteristic term, and The relationship is transposed. For variables The formula for minimizing operators, variables The solution form and formula Similar, therefore, the variable Update method and variables Same, variable The final update subproblem of is expressed as: ,in, , and the selection process of singular values ​​and variables The process of selecting singular values ​​in is the same;

[0058] When other variables remain unchanged, use simple matrix operation rules to adjust the variables To update, variables The expression of the update subproblem is:

[0059] (8);

[0060] In the formula, is the identity matrix, is the matrix multiplication rule;

[0061] Under the premise that other variables remain unchanged, the variable Update method and variables Same, both are simple matrix operations; variables The expression of the update subproblem is:

[0062] (9);

[0063] In the variable , , and Under the condition of no change, construct the solution variable The expression of the update subproblem is:

[0064] (10);

[0065] In the formula, is the Lagrangian operator of the noise term, and The relationship is transposed. For variables The formula minimization operator of ;

[0066] The soft threshold method is used to The sub-problem described is solved, and the expression of the soft threshold method is:

[0067] (11);

[0068] In the formula, is the threshold value, It is the soft threshold method;

[0069] variable The update expression is:

[0070] (12);

[0071] Mode Medium, threshold ;

[0072] variable , ,、 and The expression of the update subproblem is:

[0073]

[0074]

[0075]

[0076] (13);

[0077] In the formula, is the magnification factor, Indicates the preset The upper limit value is set to ;

[0078] The convergence conditions are: , , ;in, , , ; When the convergence condition , as well as , when both are satisfied, the iteration loop is exited, where , otherwise, continue iterating until the convergence condition is met.

[0079] Optionally, in step 3, the initial thin cloud image is subjected to feature weighting processing using the optimized thin cloud image potential low-rank representation model to obtain a low-rank optimized cloud image, including:

[0080] The basic feature items and significant feature items of the initial thin cloud map are weighted respectively to obtain a low-rank optimized cloud map;

[0081] Among them, the expression for obtaining the low-rank optimized cloud graph is:

[0082] (14);

[0083] Mode middle, is the basic feature item, is a notable feature item, Optimize the cloud map for low rank, and is the linear weight of weighted processing, and its value range is .

[0084] Optionally, in step 4, the method for calculating the building area index map and the water body index map according to the reflectivity data of multiple bands includes:

[0085] The built-up area distribution index is calculated based on the reflectance data of two short-wave infrared bands;

[0086] The normalized water index is calculated based on the reflectance data of the green band and the near-infrared band;

[0087] According to the empirical segmentation threshold, the building area distribution index and the normalized water body index are segmented to obtain the building area index map and the water body index map;

[0088] The expression of the building area distribution index is:

[0089] (15);

[0090] In the formula, is the building area distribution index of the short-wave infrared ratio, and They are respectively the reflectivity of the short-wave infrared band with a central wavelength of 1.61μm and the reflectivity of the short-wave infrared band with a central wavelength of 2.2μm;

[0091] The expression of normalized water index is:

[0092] (16);

[0093] In the formula, is the normalized water index, is the reflectivity of the green band with a central wavelength of 0.5615μm, is the reflectivity of the near-infrared band with a central wavelength of 0.8654μm;

[0094] The expressions for obtaining the building area index map and water body index map are:

[0095]

[0096] (17);

[0097] In the formula, is the threshold segmentation method, and are the empirical segmentation thresholds corresponding to SIBI and NDWI, Take the threshold The building area index map at that time, Take the threshold Water body index map at that time.

[0098] Optionally, in step 5, the method for calculating the non-water body building area distribution index map according to the building area index map and the water body index map includes:

[0099] Perform intersection operation on the building area index map and the water body index map to obtain the non-water body building area distribution index map;

[0100] Among them, the expression for obtaining the distribution index map of non-water body building areas is:

[0101] (18);

[0102] In the formula, The distribution index map of non-water body building areas Elements of This is the final distribution index map of non-water body built-up areas; 0 and 1 are the logical values ​​returned by the index, 1 is the built-up area, and 0 is the non-built-up area.

[0103] Optionally, in step 6, the method for generating a weight suppression matrix for adaptively adjusting the weight of the built-up area distribution index based on the non-water body built-up area distribution index map includes:

[0104] Based on the non-water body building area distribution index map, a non-water body building area index is obtained;

[0105] Based on the attenuation coefficient of the built-up area index and the built-up area distribution index without water interference, a weight suppression matrix for adaptively adjusting the built-up area distribution index is constructed to perform weight suppression processing on the low-rank optimization cloud map.

[0106] Among them, the expression of non-water body building area index is:

[0107] (19);

[0108] In the formula, is the non-water body built-up area index, is the index of non-water body building area, The SIBI index is the building area distribution index in the index Built-up area index at

[0109] The expression for constructing the weight suppression matrix for adaptively adjusting the built-up area distribution index is:

[0110] (20);

[0111] In the formula, is the final weight suppression matrix, is the attenuation coefficient of the building area index, ranging from [1,10], The bigger, The bigger it is, the smaller it is otherwise.

[0112] Optionally, in step 7, the low-rank optimized cloud image is weight-suppressed according to the weight suppression matrix, and the expression for obtaining the cloud image with reduced highlighted building areas is:

[0113] (twenty one);

[0114] In the formula, Optimize the cloud map for low rank, To reduce the cloud image of the highlighted building area;

[0115] It also includes the estimation of the thickness of thin clouds after lowering the cloud image of the highlighted building area.

[0116] Optionally, in step 8, the method of optimizing the cloud image with the highlighted building area reduced to obtain the final optimized thin cloud image includes:

[0117] Based on the cloud image with the highlighted building area reduced, edge search is performed to identify the edge and gradient information of thin clouds and the remaining block seams;

[0118] Nonlinear guided filtering is used to filter the edge and gradient information of thin clouds and the remaining block seams to obtain the final optimized thin cloud image, in which the block seams are smooth and the edge and gradient information of thin clouds are enhanced.

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

[0120] The present invention combines the distribution of thin clouds and building areas in remote sensing images with the low-rank sparsity theory, constructs a potential low-rank representation model for thin cloud images, and proposes a new building area index, which effectively suppresses the influence of building areas on thin cloud images. The thin cloud images can be optimized in a targeted manner for images affected by thin clouds under different time, illumination, and space factors, thereby solving the problems of poor robustness and low precision of thin cloud thickness estimation results in the prior art, and providing accurate and robust technical support for thin cloud correction and thin cloud thickness estimation of remote sensing images.

[0121] The present invention decomposes the thin cloud image into basic homogeneous cloud image features, significant non-homogeneous cloud image features and surface noise information, optimizes the thin cloud image by using adaptive suppression weights in building areas, and finally obtains the final thin cloud image through guided filtering. The final cloud image can contain less surface information, and the thin cloud intensity estimation is more accurate, which reduces the calculation amount of thin cloud thickness estimation and saves computing resources.

[0122] The present invention can also detect water bodies in the building area of ​​the image and preprocess the water body detection data into a normalized water body index to eliminate the interference of water bodies in cloudless areas on the cloud map, further reducing the impact of the building area on the estimation of thin cloud maps. BRIEF DESCRIPTION OF THE DRAWINGS

[0123] Figure 1 Shown is a flow chart of the method of the present invention;

[0124] Figure 2 Shown is a schematic diagram of thin cloud images and various bands of the present invention;

[0125] Figure 3 The figure shows the low rank analysis diagram of the thin cloud image based on SVD decomposition of the present invention;

[0126] Figure 4 Shown is a schematic diagram of different results of LatentLRR decomposition based on the present invention;

[0127] Figure 5 Shown is a distribution change diagram after building area extraction and optimization of the present invention;

[0128] Figure 6 Shown is a diagram showing the change of the thin cloud diagram optimization process of the present invention;

[0129] Figure 7 Shown is a linear relationship diagram after fitting of clean pixel samples selected from the thin cloud map of the present invention. DETAILED DESCRIPTION

[0130] The present invention will be further described below in conjunction with the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solution of the present invention, and cannot be used to limit the protection scope of the present invention.

[0131] Example 1

[0132] This embodiment provides a thin cloud thickness map optimization method considering potential low-rank features, including:

[0133] Step 1: Read a multispectral remote sensing thin cloud image, and obtain reflectance data of multiple bands in the image to estimate an initial thin cloud map; wherein the multiple bands include: a coastal aerosol band, a blue band, a green band, a red band, a near-infrared band, and two short-wave infrared bands; the coastal aerosol band ranges from 0.435 to 0.451 μm, and the blue band ranges from 0.452 to 0.512 μm; the green band ranges from 0.533 to 0.590 μm, and the red band ranges from 0.636 to 0.673 μm; the near-infrared band ranges from 0.851 to 0.879 μm; and the two short-wave infrared bands range from 1.566 to 1.651 μm and 2.107 to 2.294 μm, respectively.

[0134] Step 2: Estimate the initial thin cloud map based on the reflectivity data of multiple bands in the image;

[0135] Step 3: Use the optimized thin cloud map potential low-rank representation model to perform feature weighting processing on the initial thin cloud map to obtain a low-rank optimized cloud map.

[0136] Step 4: Based on the reflectance data of multiple bands, calculate the building area index map and water body index map;

[0137] Step 5: Calculate the distribution index map of non-water body building areas based on the building area index map and the water body index map;

[0138] Step 6: Based on the non-water body built-up area distribution index map, a weight suppression matrix for adaptively adjusting the built-up area distribution index weight is generated;

[0139] Step 7: Perform weight suppression processing on the low-rank optimized cloud image according to the weight suppression matrix to obtain a cloud image with reduced highlighted building areas;

[0140] Step 8: Optimize the cloud image that reduces the highlighted building area to obtain the final optimized thin cloud image, so as to achieve the purpose of correcting the thickness of thin clouds.

[0141] Optionally, step 2 is as follows:

[0142] According to the reflectivity data of the coastal aerosol band, the blue band, the green band, the red band, the near-infrared band and the two short-wave infrared bands, the pixels of the coastal aerosol band and the blue band in the cloudless area of ​​the image are extracted as sample data for linear fitting, and other bands are used as observation quantities for linear fitting to obtain the linear fitting relationship of the sample data in the cloudless area; in this embodiment, the central wavelength of the coastal aerosol band as the observation quantity is 0.443 μm; the central wavelength of the blue band as the observation quantity is 0.482 μm; the central wavelength of the green band as the observation quantity is 0.5615 μm; the central wavelength of the red band as the observation quantity is 0.6545 μm; the central wavelength of the near-infrared band as the observation quantity is 0.8645 μm; the central wavelengths of the two short-wave infrared bands as the observation quantities are 1.61 μm and 2.20 μm respectively;

[0143] Based on the linear fitting relationship of the sample data, the M-estimation method is used to obtain the linear relationship between the reflectance of the coastal aerosol band and the blue band. In this embodiment, two bands with central wavelengths of 0.443 μm and 0.482 μm are used;

[0144] Band synthesis is performed based on the linear relationship between the two bands with central wavelengths of 0.443μm and 0.482μm and the reflectivity of the two bands to obtain a synthetic band that is easily affected by thin clouds;

[0145] Processing the synthetic band susceptible to thin cloud influence based on a dark pixel search method, identifying dark pixels in a local window of the synthetic band image;

[0146] The dark pixels in the local window are resampled according to the bicubic interpolation method to generate an initial thin cloud image;

[0147] The expression for obtaining the synthetic band susceptible to thin clouds is:

[0148] (1);

[0149] In the formula, For the synthetic band, and They are respectively the reflectivity at the center wavelength of 0.443μm and the reflectivity at the center wavelength of 0.482μm. and is the linear relationship coefficient obtained by using the M-estimation method for the sample data in the cloud-free area.

[0150] Preferably, the optimization method of the thin cloud graph potential low-rank representation model includes:

[0151] Perform low-rank characteristic analysis on the initial thin cloud image to obtain the basic feature items and significant feature items of the initial thin cloud image; and analyze the sparse surface noise information in the initial thin cloud image;

[0152] Based on the basic feature items and significant feature items of the initial thin cloud image and the sparse surface noise information, a relationship model between low-rank thin cloud image information and sparse surface information is constructed;

[0153] According to the relationship model between low-rank thin cloud image information and sparse surface information, a convex optimization problem of features is set;

[0154] The augmented Lagrangian algorithm is used to construct the minimization optimization function for solving convex optimization problems;

[0155] Iterate multiple variables of the minimization optimization function of the convex optimization problem respectively until the convergence condition is met and stop the iteration to obtain the optimized thin cloud map potential low-rank representation model;

[0156] Among them, the expression for constructing the relationship model between low-rank thin cloud image information and sparse surface information is:

[0157] (2);

[0158] In the formula, is the initial thin cloud image, is the basic characteristic coefficient, is the basic feature item, is the significant characteristic coefficient, is a notable feature item, is the noise term, is the nuclear norm, which is used to represent The low rank property of for norm, which is used to represent Sparsity, is the regularization parameter, with a value range of [0,1];

[0159] The expression of the minimization optimization function for solving convex optimization problems is constructed as follows:

[0160] (3);

[0161] In the formula, and are the low-rank decomposition results of basic feature items and significant feature items, respectively. , and are the Lagrangian operators of noise term, basic characteristic term and significant special term, respectively. is the trace of the matrix, is the F-norm of the matrix, is the penalty parameter, .

[0162] Preferably, the method of iterating the multiple variables of the minimization optimization function of the convex optimization problem respectively until the iteration is stopped when the convergence condition is satisfied comprises:

[0163] Based on minimizing the variables of the optimization function, the optimization is solved, and the update sub-problems of each variable are established respectively. In the process of solving the update sub-problems of each variable, other variables are kept unchanged, and the update of the variable is stopped when the associated variable meets the convergence condition;

[0164] Among them, the variables include: low-rank decomposition results of basic feature items , the low-rank decomposition results of the significant feature items , basic characteristic coefficient , significant characteristic coefficient , noise term , the Lagrangian operator of the noise term , the Lagrangian operator of the basic characteristic term and the Lagrangian operator with significant special terms ;

[0165] The sub-problem solving process for each variable includes:

[0166] Keeping other variables unchanged, the variable The expression of the update subproblem is:

[0167] (4);

[0168] In the formula, For variables The formula for minimizing the operator, is the Lagrangian operator of the basic characteristic term, and They are in transposition relationship with each other;

[0169] variable The expressions of the update process are:

[0170]

[0171] (5);

[0172] In the formula, for Temporary variables, is the singular value decomposition method, are the three sub-matrices obtained after SVD decomposition, They are the left singular matrix, singular value matrix and right singular matrix respectively.

[0173] Using the singular value threshold method, the singular values ​​after SVD decomposition are Screen and get the number of singular values ​​that meet the conditions , in this way Update and get the expression of the number of singular values ​​that meet the conditions:

[0174]

[0175] (6);

[0176] In the formula, is the number of singular values ​​that meet the conditions, is a singular value vector;

[0177] Get variables The final update subproblem of is expressed as:

[0178]

[0179] In the formula, is a diagonal matrix;

[0180] Keeping other variables unchanged, the variable The expression of the update subproblem is:

[0181] (7);

[0182] In the formula, is the Lagrangian operator of the significant characteristic term, and The relationship is transposed. For variables The formula for minimizing operators, variables The solution form and formula Similar, therefore, the variable Update method and Same, variable The final update subproblem of is expressed as: ,in, , and the selection process of singular values ​​and variables The process of selecting singular values ​​in is the same;

[0183] When other variables remain unchanged, use simple matrix operation rules to adjust the variables To update, variables The expression of the update subproblem is:

[0184] (8);

[0185] In the formula, is the identity matrix, is the matrix multiplication rule;

[0186] Under the premise that other variables remain unchanged, the variable Update method and variables Same, both are simple matrix operations; variables The expression of the update subproblem is:

[0187] (9);

[0188] In the variable , , and Under the condition of no change, construct the solution variable The expression of the update subproblem is:

[0189] (10);

[0190] In the formula, is the Lagrangian operator of the noise term, and The relationship is transposed. For variables The formula minimization operator of ;

[0191] The soft threshold method is used to The sub-problem described is solved, and the expression of the soft threshold method is:

[0192] (11);

[0193] In the formula, is the threshold value, It is the soft threshold method;

[0194] variable The update expression is:

[0195] (12);

[0196] Mode Medium, threshold .

[0197] variable , ,、 and The expression of the update subproblem is:

[0198]

[0199]

[0200]

[0201] (13);

[0202] In the formula, is the magnification factor, Indicates the preset The upper limit is set here as ;

[0203] The convergence conditions are: , , ;in, , , ; When the convergence condition , as well as , when both are satisfied, the iteration loop is exited, where , otherwise, continue iterating until the convergence condition is met.

[0204] Preferably, considering that there are overlapping cloud areas in the basic features and the significant features, direct addition will cause the thickness of the overlapping area to increase, so in step 3, the initial thin cloud image is subjected to feature weighting processing using the optimized thin cloud image potential low-rank representation model to obtain a low-rank optimized cloud image;

[0205] Among them, the expression for obtaining the low-rank optimized cloud graph is:

[0206] (14);

[0207] Mode middle, Optimize the cloud map for low rank, and is the linear weight of weighted processing, and its value range is .

[0208] Preferably, in step 4, the building area index map and the water body index map are calculated based on the reflectivity data of multiple bands, including:

[0209] The built-up area distribution index is calculated based on the reflectance data of two short-wave infrared bands;

[0210] The normalized water index is calculated based on the reflectance data of the green band and the near-infrared band;

[0211] According to the empirical segmentation threshold, the building area distribution index and the normalized water body index are segmented to obtain the building area index map and the water body index map;

[0212] The expression of the building area distribution index is:

[0213] (15);

[0214] In the formula, is the building area distribution index of the short-wave infrared ratio, and They are respectively the reflectivity of the short-wave infrared band with a central wavelength of 1.61μm and the reflectivity of the short-wave infrared band with a central wavelength of 2.2μm;

[0215] The expression of normalized water index is:

[0216] (16);

[0217] In the formula, is the normalized water index, is the reflectivity of the green band with a central wavelength of 0.5615μm, is the reflectivity of the near-infrared band with a central wavelength of 0.8654μm;

[0218] The expressions for obtaining the building area index map and water body index map are:

[0219]

[0220] (17);

[0221] In the formula, is the threshold segmentation method, and are the empirical segmentation thresholds corresponding to SIBI and NDWI, Take the threshold The building area index map at that time, Take the threshold Water body index map at that time.

[0222] Step 5 includes:

[0223] Perform intersection operation on the building area index map and the water body index map to obtain the non-water body building area distribution index map;

[0224] Among them, the expression for obtaining the distribution index map of non-water body building areas is:

[0225] (18);

[0226] In formula (18), The distribution index map of non-water body building areas Elements of This is the final distribution index map of non-water body built-up areas; 0 and 1 are the logical values ​​returned by the index, 1 is the built-up area, and 0 is the non-built-up area.

[0227] Preferably, step 6 comprises: obtaining a non-water body building area index based on the non-water body building area distribution index map;

[0228] Based on the attenuation coefficient of the built-up area index and the built-up area distribution index without water interference, a weight suppression matrix for adaptively adjusting the built-up area distribution index is constructed to perform weight suppression processing on the low-rank optimization cloud map.

[0229] Among them, the expression of non-water body building area index is:

[0230] (19);

[0231] In the formula, is the non-water body built-up area index, is the index of non-water body building area, The SIBI index is the building area distribution index in the index Built-up area index at

[0232] The expression for constructing the weight suppression matrix for adaptively adjusting the built-up area distribution index is:

[0233] (20);

[0234] In the formula, is the final weight suppression matrix, and ; is the attenuation coefficient of the built-up area index (SIBI), ranging from [1,10], The bigger, The bigger it is, the smaller it is otherwise.

[0235] Preferably, in step 7, the low-rank optimized cloud image is weighted suppressed according to the weight suppression matrix, and the expression for obtaining the cloud image with reduced highlighted building areas is:

[0236] (twenty one);

[0237] In the formula, To reduce the cloud image of the highlighted building area;

[0238] In addition, step 7 also includes obtaining a thin cloud thickness estimation result after lowering the cloud image of the highlighted building area. When lowering the cloud image of the highlighted building area, it is possible to obtain estimation information of the thickness of all thin clouds and other parameters of thin clouds in the cloud area, thereby predicting and evaluating the thin cloud thickness data without interference from the building area, water body and other factors.

[0239] Preferably, step 8 includes: performing edge search based on the cloud image with the highlighted building area reduced, identifying the edge and gradient information of the thin cloud and the remaining block seams;

[0240] Nonlinear guided filtering is used to filter the edge and gradient information of thin clouds and the residual block seams, thereby smoothing the block seams and enhancing the edge and gradient information of thin clouds to obtain the final optimized thin cloud image.

[0241] Example 2

[0242] This embodiment provides a thin cloud thickness map optimization method considering potential low-rank features, combined with Figure 1 Shown include:

[0243] Step 1: Read the multispectral remote sensing thin cloud image to obtain the reflectivity data of multiple bands in the image.

[0244] Step 2: Estimating the initial thin cloud map, constructing the building area distribution index, and constructing the normalized water body index respectively; performing feature weighting processing on the initial thin cloud map based on the optimized potential low-rank representation model of the thin cloud map to obtain a low-rank optimized cloud map;

[0245] The building area distribution index is obtained by substituting the reflectance data of the short-wave infrared bands with central wavelengths of 1.61 μm and 2.2 μm into the building area distribution index formula (15);

[0246] The normalized water index is based on the central wavelength. and Substitute the reflectance data of the green band and near-infrared band into the normalized water index Formula (16) is used to obtain:

[0247] According to the building area distribution index and the normalized water body index, the building area index map and the water body index map are calculated.

[0248] Step 3: Calculate the distribution index map of non-water body building areas based on the building area index map and the water body index map;

[0249] Step 4: Based on the non-water body built-up area distribution index map, a weight suppression matrix for adaptively adjusting the built-up area distribution index weight is generated;

[0250] Step 5: Perform weight suppression processing on the low-rank optimized cloud image according to the weight suppression matrix to obtain a cloud image with reduced highlighted building areas;

[0251] Step 6: Optimize the cloud image that reduces the highlighted building area to obtain the final optimized thin cloud image.

[0252] Combination Figure 2 and Figure 7 Step 1 of this embodiment is further described as follows:

[0253] Step 1.1: Select clean pixels in the cloud-free area of ​​the image as sample data for the linear fitting relationship, and use the M-estimation method to calculate the linear weights of the two bands with central wavelengths of 0.443μm and 0.482μm in the thin cloud image. For the selection of samples, pixels that are visually cloudless in the thin cloud image and cannot be pixels of a single type can be selected, such as Figure 2 The red area in (a) shows the pixels selected in this example.

[0254] Step 1.2: The linear relationship after fitting the clean pixel samples selected in step 1.1 is as follows. In this example, the linear weight coefficients solved based on the selected samples are a=0.7571, b=0.0461, and the fitting accuracy R2=0.9904.

[0255] When clouds appear in the scene, the fitted line of this clean pixel will have a large deviation, based on which a new band susceptible to the influence of thin clouds is synthesized;

[0256] In this example, the newly synthesized bands are Figure 2 (d) The synthetic band map shown.

[0257] Step 1.3: According to the dark pixel search method, such as non-overlapping minimum filtering, the synthetic band obtained by formula (1) is Perform a dark pixel search with the window set to 3, and then use the bicubic interpolation method to resample the result to the original size to obtain the initial thin cloud image, as shown in Figure 3 (a).

[0258] The estimation of the initial thin cloud image in step 2 also includes:

[0259] Step 2.1: Construct a relationship model between low-rank thin cloud image information and sparse surface information. For the initial thin cloud image, it has the characteristics of low rank. Figure 3 (a) to (g) shown in the figure are all based on SVD decomposition. Figure 3 (a) is the initial thin cloud image, Figure 3 (b) is the left singular matrix U diagram, Figure 3 (c) is the right singular matrix V T picture, Figure 3 (d) is the singular value matrix S graph, Figure 3 (e) is N S = 50, Figure 3 (f) is N S = 30, Figure 3 (g) is N S = 10; where N S is the number of singular values ​​on the diagonal of the singular value matrix S.

[0260] After low-rank analysis of the initial thin cloud image, basic feature items and significant feature items are obtained, namely, homogeneous thin clouds and prominent thin clouds. In addition, it also contains sparse surface noise information, such as building areas.

[0261] Step 2.2: Use the augmented Lagrangian algorithm to construct a minimization optimization model. The characteristics of the model can be identified as a convex optimization problem. For the optimization solution of convex optimization problems, the classic augmented Lagrangian algorithm can be used.

[0262] Step 2.3: Optimize and solve the constructed minimization optimization model, specifically: decompose the variables involved in formula (3) into multiple sub-problems for solution. In this embodiment, multiple variables are decomposed into 6 sub-problems, and the sub-problem expressions are the same as formulas (4) to (13) in Example 1;

[0263] The convergence conditions are: , , ;in, , , ; When the convergence condition , as well as , when both are satisfied, the iteration loop is exited, where , otherwise, continue iterating until the convergence condition is met.

[0264] The basic feature items and significant feature items of the initial thin cloud image contain non-homogeneous thin cloud features. Based on this, in order to avoid over-correction, a linear weighted method is used to obtain a low-rank optimized thin cloud image, which can be expressed as . Figure 4 Shown when , , , the decomposed basic feature items, significant feature items, and noise items are shown in Figure 4 (a)–(c). In addition, the linear weighting coefficient is set and Get the final low-rank optimized cloud map, see Figure 4 (d).

[0265] Step 2.5: Use the reflectance data of the short-wave infrared bands with central wavelengths of 1.61μm and 2.2μm to construct the building area distribution index associated with the cloud map. Due to the influence of thin clouds, the extraction of the building area will be greatly affected. As the wavelength increases, its penetration becomes stronger and the influence of thin clouds becomes weaker. Based on this, a building area distribution index based on the short-wave infrared ratio (Short-wave Infrared based Building Index, referred to as SIBI) is proposed, which can be expressed as formula (15), as follows Figure 5 (a) shows the index diagram of the built-up area of ​​this embodiment, from which it can be seen that the value of the built-up area is significantly higher than that of other areas.

[0266] Step 3: Get the distribution index map of non-water body building areas , since the influence of the built-up area on the estimation of thin cloud images cannot be ignored, after obtaining the built-up area distribution index through step 2, the built-up area and non-built-up area can be segmented according to the threshold. In order to prevent the built-up area from containing water bodies, the Normalized Difference Water Index (NDWI) is used to reprocess it to eliminate the interference of water bodies in the cloudless area, which can be expressed as formula (16), as follows Figure 5 (c) shows the NDWI diagram. In formula (17) and formula (18), and are the empirical segmentation thresholds corresponding to SIBI and NDWI respectively. Here, for The built-up area index map at the time of , which consists of logical values ​​0 and 1, 1 represents the built-up area, and 0 represents the non-built-up area, such as Figure 5 (b) As shown. Similarly, for Water body index map at time Then it is the final non-water building area distribution index map, such as Figure 5 (d) shown.

[0267] In step 4, based on the non-water body building area distribution index map, a weight suppression matrix for adaptively adjusting the weight of the building area distribution index is generated to re-optimize the low-rank optimized cloud map obtained in step 2, thereby suppressing the remaining highlighted building areas in the cloud map and further reducing the impact of the building areas on the low-rank optimized cloud map. The calculation process can be expressed by equations (19)-(21).

[0268] In this example, when the attenuation coefficient of the built-up area index (SIBI) When , the requirement can be met. The minimum value of is 0.78, so The maximum value of , which is less than 1, thus achieving the suppression effect on the thin cloud map in the building area. The results are shown in Figure 6 (c). Figure 6 (a) is the initial thin cloud image, Figure 6 (b) is the low-rank optimized cloud map, Figure 6 (c) is a cloud image to reduce the highlight of the building area; Figure 6 (d) is the final optimized thin cloud image. Figure 6 In (a)-(d), the adjustment effect can be clearly seen. If it is used in other scenarios, it can be adjusted according to needs.

[0269] Step 6 obtains the final optimized thin cloud image, wherein the block seam effect remaining in the building area in the final optimized thin cloud image is smoothed by using nonlinear guided filtering, and the edge and gradient information of the thin cloud image is retained, thereby obtaining smooth block seams and enhanced thin cloud edge and gradient information, thereby achieving the purpose of clarity correction of the thin cloud image and thin cloud thickness estimation, as shown in FIG. Figure 6 (d) shows an example of the final optimized thin cloud map.

[0270] In summary, the present invention combines the distribution of thin clouds and building areas in remote sensing images with the low-rank sparsity theory, constructs a potential low-rank representation model of thin cloud images, and proposes a new type of building area index, which effectively suppresses the influence of building areas on thin cloud images. It can optimize thin cloud images in a targeted manner for images affected by thin clouds under different time, illumination, and space factors, solve the problems of poor robustness and low precision of thin cloud thickness estimation results in the prior art, and provide accurate and robust technical support for thin cloud correction and thin cloud thickness estimation of remote sensing images. Secondly, the present invention decomposes the thin cloud image into basic homogeneous cloud image features, significant non-homogeneous cloud image features, and surface noise information, and optimizes the thin cloud image using adaptive suppression weights of building areas, and then obtains the final thin cloud image through guided filtering. The final cloud image can contain less surface information, and the thin cloud intensity estimation is more accurate, which reduces the amount of calculation for thin cloud thickness estimation and saves computing resources.

[0271] The present invention can also detect water bodies in the building area of ​​the image and preprocess the water body detection data into a normalized water body index to eliminate the interference of water bodies in cloudless areas on the cloud map, further reducing the impact of the building area on the estimation of thin cloud maps.

[0272] It will be appreciated by those skilled in the art that embodiments of the present invention may be provided as methods, systems, or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Furthermore, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0273] The present invention is described with reference to flowcharts and / or block diagrams of methods, devices (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0274] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture including an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0275] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process in the computer or other programmable device. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0276] The embodiments of the present invention are described above in conjunction with the accompanying drawings, but the present invention is not limited to the above-mentioned specific implementation methods. The above-mentioned specific implementation methods are merely illustrative and not restrictive. Under the enlightenment of the present invention, ordinary technicians in this field can also make many forms without departing from the scope of protection of the purpose of the present invention and the claims, which all fall within the protection of the present invention.

Claims

1. A method for optimizing thin cloud thickness maps considering potential low-rank features, characterized in that: include: Read multispectral remote sensing thin cloud images; An initial thin cloud map is estimated based on the reflectivity data of multiple bands in the image; the multiple bands include: coastal aerosol band, blue band, green band, red band, near infrared band and two short-wave infrared bands; Using an optimized potential low-rank representation model of thin cloud images, feature weighting is performed on the initial thin cloud image to obtain a low-rank optimized cloud image; Based on the reflectance data of multiple bands, the building area index map and water body index map are calculated; According to the building area index map and the water body index map, the distribution index map of the non-water body building area is calculated; Based on the non-water body built-up area distribution index map, a weight suppression matrix for adaptively adjusting the built-up area distribution index weight is generated; The low-rank optimized cloud image is weight-suppressed according to the weight suppression matrix to obtain a cloud image with reduced highlighted building areas; Optimize the cloud image that reduces the highlighted building area to obtain the final optimized thin cloud image; The optimization methods of the thin cloud graph potential low-rank representation model include: Perform low-rank characteristic analysis on the initial thin cloud image to obtain the basic feature items and significant feature items of the initial thin cloud image; and analyze the sparse surface noise information in the initial thin cloud image; Based on the basic feature items and significant feature items of the initial thin cloud image and the sparse surface noise information, a relationship model between low-rank thin cloud image information and sparse surface information is constructed; According to the relationship model between low-rank thin cloud image information and sparse surface information, a convex optimization problem of features is set; The augmented Lagrangian operator method is used to construct the minimization optimization function for solving convex optimization problems; Iterate multiple variables of the minimization optimization function of the convex optimization problem respectively until the convergence condition is met and stop the iteration to obtain the optimized thin cloud map potential low-rank representation model; Among them, the expression of the relationship model between low-rank thin cloud image information and sparse surface information is: （1）； In the formula, is the initial thin cloud image, is the basic characteristic coefficient, is the basic feature item, is the significant characteristic coefficient, is a notable feature item, is the noise term, is the nuclear norm, indicating low rank, for norm, indicating sparsity, is the regularization parameter, with a value range of [0,1]; The expression of the minimized optimization function for solving the convex optimization problem is: （2）； In the formula, and are the low-rank decomposition results of basic feature items and significant feature items, respectively. , and are the Lagrangian operators of noise term, basic characteristic term and significant special term, respectively. is the trace of the matrix, is the F-norm of the matrix, is the penalty parameter, ; The initial thin cloud image is subjected to feature weighting processing using the optimized thin cloud image potential low-rank representation model to obtain a low-rank optimized cloud image, including: The basic feature items and significant feature items of the initial thin cloud map are weighted respectively to obtain a low-rank optimized cloud map; Among them, the expression for obtaining the low-rank optimized cloud graph is: （3）； In the formula, is the basic feature item, is a notable feature item, Optimize the cloud map for low rank, and They are all linear weights for weighted processing, and their value range is .

2. The thin cloud thickness map optimization method considering potential low-rank features according to claim 1 is characterized in that: The initial thin cloud map is estimated based on the reflectivity data of multiple bands in the image, including: Obtain the reflectance data of the coastal aerosol band, blue band, green band, red band, near infrared band and two short-wave infrared bands, and extract the pixels of the coastal aerosol band and blue band in the cloud-free area of ​​the image as the sample data of the linear fitting, and use other bands as the observed quantities of the linear fitting to obtain the linear fitting relationship of the sample data in the cloud-free area; Based on the linear fitting relationship of the sample data, the M-estimation method is used to obtain the linear relationship between the reflectance of the coastal aerosol band and the blue band; Band synthesis is performed based on the linear relationship and reflectivity between the coastal aerosol band and the blue band to obtain a synthetic band that is easily affected by thin clouds; Processing the synthetic band susceptible to thin cloud influence based on a dark pixel search method, identifying dark pixels in a local window of the synthetic band image; The dark pixels in the local window are resampled according to the bicubic interpolation method to generate an initial thin cloud image; The expression for obtaining the synthetic band susceptible to thin clouds is: （4）； In the formula, For the synthetic band, and are the reflectances of the coastal aerosol band and the blue band, and These are the linear relationship coefficients obtained using the M-estimation method for sample data in the cloud-free area.

3. The thin cloud thickness map optimization method considering potential low-rank features according to claim 1 is characterized in that: Iterate multiple variables of the minimization optimization function of the convex optimization problem respectively until the iteration is stopped when the convergence condition is met, including: Based on minimizing the variables of the optimization function, the optimization is solved, and the update sub-problems of each variable are established respectively. In the process of solving the update sub-problems of each variable, other variables are kept unchanged, and the update of the variable is stopped when the associated variable meets the convergence condition; Among them, the variables include: low-rank decomposition results of basic feature items , the low-rank decomposition results of the significant feature items , basic characteristic coefficient , significant characteristic coefficient , noise term , the Lagrangian operator of the noise term , the Lagrangian operator of the basic characteristic term and the Lagrangian operator with significant special terms ; The sub-problem solving process for each variable includes: Keeping other variables unchanged, the variable The expression of the update subproblem is: （5）； In the formula, For variables The formula for minimizing the operator, is the Lagrangian operator of the basic characteristic term, and They are in transposition relationship with each other; variable The expressions of the update process are: （6）； In the formula, for Temporary variables, is the singular value decomposition method, are the three sub-matrices obtained after SVD decomposition, They are the left singular matrix, singular value matrix and right singular matrix respectively; Using the singular value threshold method, the singular values ​​after SVD decomposition are Screen and get the number of singular values ​​that meet the conditions , in this way Update and get the expression of the number of singular values ​​that meet the conditions: （7）； In the formula, is the number of singular values ​​that meet the conditions, is a singular value vector; Get variables The final update subproblem of is expressed as: In the formula, is a diagonal matrix; Keeping other variables unchanged, the variable The expression of the update subproblem is: （8）； In the formula, is the Lagrangian operator of the significant characteristic term, and The relationship is transposed. For variables The formula for minimizing operators, variables Update method and same; When other variables remain unchanged, To update, variables The expression of the update subproblem is: （9）； In the formula, is the identity matrix, is the matrix multiplication rule; Under the premise that other variables remain unchanged, To update, variables The expression of the update subproblem is: （10）； In the variable , , and Under the condition of no change, construct the solution variable The update subproblem expression is: （11）； In the formula, is the Lagrangian operator of the noise term, and The relationship is transposed. For variables The formula minimization operator of ; The soft threshold method is used to The sub-problem described is solved, and the expression of the soft threshold method is: （12）； In the formula, is the threshold value, It is the soft threshold method; variable The update expression is: （13）； In the formula, the threshold ; variable , ,、 and The expression of the update subproblem is: （14）； In the formula, is the magnification factor, For pre-set The upper limit is set to ; The convergence conditions are: , , ;in, , , ; When the convergence conditions are met at the same time , as well as When , it jumps out of the iteration loop, where , otherwise, continue iterating until the convergence condition is met.

4. The thin cloud thickness map optimization method considering potential low-rank features according to claim 1 is characterized in that: Based on the reflectance data of multiple bands, the building area index map and water body index map are calculated, including: The built-up area distribution index is calculated based on the reflectance data of two short-wave infrared bands; The normalized water index is calculated based on the reflectance data of the green band and the near-infrared band; According to the empirical segmentation threshold, the building area distribution index and the normalized water body index are segmented to obtain the building area index map and the water body index map; The expression of the building area distribution index is: （15）； In the formula, is the building area distribution index of the short-wave infrared ratio, and They are respectively the reflectivity of the short-wave infrared band with a central wavelength of 1.61μm and the reflectivity of the short-wave infrared band with a central wavelength of 2.2μm; The expression of normalized water index is: （16）； In the formula, is the normalized water index, is the reflectivity of the green band with a central wavelength of 0.5615μm, is the reflectivity of the near-infrared band with a central wavelength of 0.8654μm; The expressions for obtaining the building area index map and water body index map are: （17）； In the formula, is the threshold segmentation method, and are the empirical segmentation thresholds corresponding to SIBI and NDWI, Take the threshold The building area index map at that time, Take the threshold Water body index map at that time.

5. The thin cloud thickness map optimization method considering potential low-rank features according to claim 4 is characterized in that: According to the building area index map and water body index map, the non-water body building area distribution index map is calculated, including: Perform intersection operation on the building area index map and the water body index map to obtain the non-water body building area distribution index map; Among them, the expression for obtaining the distribution index map of non-water body building areas is: （18）； In the formula, The distribution index map of non-water body building areas Elements of This is the final distribution index map of non-water body built-up areas; 0 and 1 are the logical values ​​returned by the index, 1 is the built-up area, and 0 is the non-built-up area.

6. The thin cloud thickness map optimization method considering potential low-rank features according to claim 5 is characterized in that: Based on the non-water body built-up area distribution index map, a weight suppression matrix for adaptively adjusting the built-up area distribution index weight is generated, including: Based on the non-water body building area distribution index map, a non-water body building area index is obtained; Based on the attenuation coefficient of the built-up area index and the non-water body built-up area index, a weight suppression matrix for adaptively adjusting the built-up area distribution index is constructed to perform weight suppression processing on the low-rank optimization cloud map. Among them, the expression of non-water body building area index is: （19）； In the formula, is the non-water body built-up area index, is the index of non-water body building area, The SIBI index is the building area distribution index in the index The built-up area index at The expression for constructing the weight suppression matrix for adaptively adjusting the built-up area distribution index is: （20）； In the formula, is the final weight suppression matrix, is the attenuation coefficient of the building area index, and its value range is [1,10].

7. The thin cloud thickness map optimization method considering potential low-rank features according to claim 6 is characterized in that: The low-rank optimized cloud image is weighted suppressed according to the weight suppression matrix, and the expression of the cloud image with reduced highlighted building areas is obtained as follows: （21）； In the formula, Optimize the cloud map for low rank, This is the cloud map result after the weight of the building area is optimized; It also includes the estimation of the thickness of thin clouds after lowering the cloud image of the highlighted building area.

8. The method for optimizing the thin cloud thickness map considering potential low-rank features according to claim 1, characterized in that: The cloud image with reduced highlight building areas is optimized to obtain the final optimized thin cloud image, including: Based on the cloud image with the highlighted building area reduced, edge search is performed to identify the edge and gradient information of thin clouds and the remaining block seams; Nonlinear guided filtering is used to filter the edge and gradient information of thin clouds and the remaining block seams to obtain the final optimized thin cloud image, in which the block seams are smoothed and the edge and gradient information of thin clouds are enhanced.

Images