A remote sensing image cloud area reconstruction method based on SAR prior knowledge guidance

Abstract

The application discloses a remote sensing image cloud area reconstruction method based on SAR prior knowledge guidance. The method firstly inputs the obtained optical remote sensing image into a cloud detection module to obtain a preliminary cloud detection result, then constructs a SAR feature extraction module, inputs the obtained adjacent time phase SAR image into the SAR feature extraction module, constructs a polarization covariance matrix of the SAR image, extracts relevant polarization scattering features, finally constructs a cloud area image reconstruction module, and sequentially inputs the obtained optical remote sensing image, the cloud detection result and the polarization scattering features into the cloud area image reconstruction module, so that the spatial structure feature reconstruction and the global consistency feature repair of the cloud-shielded area of the optical remote sensing image are realized, and finally a high-fidelity cloud-free image is generated. The application effectively improves the feature recovery capability of the cloud-shielded area of the optical remote sensing image by fusing the SAR image texture and polarization features, and solves the problems of optical remote sensing image information blurring and partial information loss caused by serious cloud and fog interference.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of remote sensing image processing, and particularly relates to a remote sensing image cloud area reconstruction method based on SAR prior knowledge guidance. BACKGROUND

[0002] Optical remote sensing images have rich spectral information and can intuitively reflect the characteristics of different types of ground objects. However, due to the imaging mechanism, the ground information is easily lost due to cloud cover, which seriously affects the interpretation and interpretation of the image. According to relevant research estimates, the cloud coverage in the atmosphere is 67%, which means that optical remote sensing will inevitably be contaminated by clouds during imaging. Whether it is long-term crop phenology observation or short-term natural disaster warning, cloud-contaminated images will affect the accuracy and timeliness of monitoring. Therefore, it is necessary to design an optical remote sensing image cloud area reconstruction method to improve the usability of the image.

[0003] Synthetic aperture radar (SAR) is an active remote sensing system with a microwave band working frequency, which can actively emit signals to detect land targets, and thus has all-weather and all-day working ability, which can effectively avoid cloud pollution. Therefore, the method of using SAR images to reconstruct the cloud area of optical remote sensing images has attracted widespread attention. The existing methods of fusing SAR images to remove clouds from optical images can be divided into two categories: based on prior information and based on data generation. The method based on prior information mainly extracts the texture, geometry and other spatial detail information of the SAR image as prior constraint conditions to assist the reconstruction of the cloud-shielded part of the optical image. The data generation method directly uses the generative adversarial network (GAN) to convert the SAR image into an optical image, which is then supplemented to the cloud missing part of the optical image. Both methods have achieved good repair results, but the reconstructed image still has problems such as blurring, artifacts, rough texture, and other problems, making it difficult to perfectly reconstruct the missing ground details under cloud cover. SUMMARY

[0004] In view of the deficiencies of the prior art, the application discloses a remote sensing image cloud area reconstruction method based on SAR prior knowledge guidance, which can comprehensively consider the imaging differences between heterogeneous remote sensing images, use the prior information provided by the SAR image, and construct an optical remote sensing image cloud-shielded area repair network based on SAR prior information assistance to realize high fidelity of the cloud area reconstruction image in color, texture, and details, and finally obtain a cloud-shielded image.

[0005] In order to achieve the above purpose, the application provides a remote sensing image cloud area reconstruction method based on SAR prior knowledge guidance, comprising the following steps:

[0006] Step 1, constructing a cloud detection module, inputting the acquired optical remote sensing image into the cloud detection module for cloud spectrum testing and cloud probability calculation to obtain a cloud detection result;

[0007] Step 1.1, performing cloud spectrum testing on the acquired optical remote sensing image to obtain a preliminary cloud detection result;

[0008] Step 1.2, extracting a cloud edge zone in the image by calculating a cloud probability to obtain a cloud edge zone detection result;

[0009] Step 1.3, fusing the preliminary cloud detection result and the cloud edge zone detection result to obtain a final optical remote sensing image cloud detection result;

[0010] Step 2, constructing a SAR feature extraction module, inputting the acquired near-time SAR image into the SAR feature extraction module, constructing a polarized covariance matrix of the SAR image, and extracting relevant polarization scattering features;

[0011] Step 2.1, constructing a polarized covariance matrix of the SAR image;

[0012] Step 2.2, decomposing the polarized covariance matrix into a sum of relevant projections and irrelevant projections in different directions;

[0013] Step 2.3, calculating the polarization scattering entropy, average scattering angle, and polarization scattering inverse entropy using the polarized covariance matrix decomposition result;

[0014] Step 3, constructing a cloud area image reconstruction module, inputting the acquired optical remote sensing image, the cloud detection result obtained in step 1, and the polarization scattering features obtained in step 2 into the cloud area image reconstruction module in sequence to realize spatial structure feature reconstruction and global consistency feature repair of the cloud-shielded area of the optical remote sensing image, and finally generate a high-fidelity cloud-free image;

[0015] Step 3.1, using a local fine texture feature encoding submodule to map the polarization scattering entropy, average scattering angle, and polarization scattering inverse entropy extracted in step 2 to a feature space;

[0016] Step 3.2, sequentially outputting the last four features from the local fine texture feature encoding submodule last intermediate layer and three decoder layers S k ,k∈k0,1,2,3} from coarse to fine;

[0017] Step 3.3, based on the zero initialization residual addition feature fusion technology, gradually adding the four features S k ,k∈k0,1,2,3} obtained in step 3.2 to the global attention & receptive field joint enhancement repair submodule from coarse to fine;

[0018] Step 3.4, inputting the cloud optical remote sensing image and the cloud detection result obtained in step 1 into a global attention & receptive field joint enhancement repair sub-module to realize repair of the cloud-shielded area and other areas not shielded by the cloud in terms of global consistency, and obtain a final cloud-removed image.

[0019] Moreover, the step 1.1 specifically comprises the following sub-steps:

[0020] Step 1.1.1, performing basic testing on the optical remote sensing image;

[0021] Let the preliminary cloud detection result be R a , R a is a binary image, wherein 1 represents a cloud pixel, and 0 represents a non-cloud pixel, the short-wave infrared 2 band threshold is set as γ1, the normalized vegetation index threshold is set as γ2, and the normalized snow index threshold is set as γ3, to obtain the preliminary cloud detection result R a1 of the optical remote sensing image, that is:

[0022]

[0023]

[0024]

[0025] In the formula, ρ Red , ρ Green , and ρ SWIR2 respectively represent the reflectivity of the red light band, the green light band, and the short-wave infrared 2 band, NDVI is the normalized vegetation index, and NDSI is the normalized snow index.

[0026] Step 1.1.2, performing whiteness testing on the optical remote sensing image;

[0027] The whiteness testing is performed by calculating the difference between each visible light band and the average value of the visible light band divided by the average value of the visible light band, that is:

[0028]

[0029] meanVis=(ρ Red +ρ Green +ρ Blue ) / 3 (5)

[0030] In the formula, whiteness represents the whiteness testing, meanVis represents the average value of each visible light band, and ρ Red , ρ Green , and ρ Blue are respectively the reflectivity of the red band, the green band, and the blue light band.

[0031] The threshold of the whiteness test result is set as γ4, and a preliminary cloud detection result R is obtained a2 = whiteness < γ4

[0032] Step 1.1.3, haze detection is performed on the optical remote sensing image;

[0033] The haze optimization conversion HOT algorithm is used to set the score greater than 0 as a cloud pixel, and a preliminary cloud detection result R is obtained a3 = HOT > 0, and the specific calculation formula of HOT is as follows:

[0034] HOT = p Blue -0.5 x p Red -0.08 (6)

[0035] In the formula, p Red , p Blue are the reflectivities of the red band and the blue band respectively;

[0036] Step 1.1.4, rock and desert detection is performed on the optical remote sensing image;

[0037] By calculating the ratio of the near-infrared band and the short-wave infrared 1 band of the image, the pixels with a result greater than γ5 are set as cloud pixels, and a preliminary cloud detection result R is obtained a4 , that is:

[0038] R a4 = p NIR / p SWIR1 > γ5 (7)

[0039] In the formula, p NIR , p SWIR1 are the reflectivities of the near-infrared band and the short-wave infrared 1 band respectively;

[0040] Step 1.1.5, cirrus cloud test is performed on the optical remote sensing image;

[0041] The reflectivity of the cirrus cloud band is divided by 0.04, and the pixels with a result greater than γ6 are set as cloud pixels, and a preliminary cloud detection result R is obtained a5 , that is

[0042] R a5 = p Cirrus / 0.04 > γ6 (8)

[0043] In the formula, p Cirrus represents the reflectivity of the cirrus cloud band;

[0044] Step 1.1.6, the five preliminary cloud detection results obtained in steps 1.1.1-1.1.5 are fused to obtain a final preliminary cloud detection result;

[0045] Final preliminary cloud detection results R a Represented as:

[0046] R a =R a1 ∪R a2 ∪R a3 ∪R a4 ∪R a5 (9)

[0047] Step 1.1.7: Calculate the difference in spectral reflectance between cloud and land pixels to obtain the probability of spectral variation;

[0048] Spectral change probability L vari The calculation formula is as follows:

[0049] L vari =1-max(|NDVI|, |NDSI|, |whitenessD (10)

[0050] In the formula, NDVI is the normalized vegetation index, NDSI is the normalized snow cover index, and whiteness is the whiteness test result.

[0051] Step 1.1.8: Normalize the HOT results obtained in Step 1.1.3;

[0052] The normalization calculation formula is as follows:

[0053]

[0054] In the formula, L HOT This represents the normalized HOT result. low and HOT high These represent the 17.5th and 82.5th percentiles of the HOT value, respectively.

[0055] Furthermore, step 1.2 involves combining the preliminary cloud detection results obtained in step 1.1.5, the spectral change probability index obtained in step 1.1.7, and the HOT normalization result obtained in step 1.1.8 to calculate the cloud probability. p Cloud p Values ​​greater than γ7 are set as cloud pixels, and the cloud edge detection result R is obtained. b Cloud edge detection results R b It is a binary image where 1 represents a cloud pixel, 0 represents a non-cloud pixel, and R b The calculation method is as follows:

[0056] R b =Cloud p =L vari ×LHoT +R a5 ×0.5>γ7 (12)

[0057] where R a5 is the preliminary cloud detection result, L vari is the probability of spectral variation, L HOT is the HOT normalized result.

[0058] Moreover, the step 1.3 is to fuse the preliminary cloud detection result obtained in step 1.1.5 and the cloud edge zone detection result obtained in step 1.2, to obtain the final cloud detection result R cloud of the optical remote sensing image, where R cloud is a binary image, 1 represents a cloud pixel, and 0 represents a non-cloud pixel. cloud The calculation method is as follows:

[0059] R cloud = R a ∩ R b (13)

[0060] where R a is the preliminary cloud detection result, and R b is the cloud edge zone detection result.

[0061] Moreover, the polarization covariance matrix C2 in step 2.1 is constructed as follows:

[0062]

[0063] where V represents the vertical polarization state of the SAR antenna; H represents the horizontal polarization state of the SAR antenna; S VH represents the relationship between the incident field and the scattering field in the polarization electromagnetic wave state of transmitting a vertical signal and receiving a horizontal signal; S VV represents the relationship between the incident field and the scattering field in the polarization electromagnetic wave state of transmitting a vertical signal and receiving a vertical signal; * represents complex conjugate; <·> represents statistical average; |S VH | and |S VV | represent the phase information of the SAR image in different polarization states, and i is an imaginary unit; represents the polarization information of the SAR image.

[0064] Moreover, the step 2.2 of the eigenvalue decomposition of the polarization covariance matrix C2 obtains the eigenvectors and eigenvalues, and the eigenvectors are sorted according to the corresponding eigenvalue size to obtain a set of principal components, wherein the first principal component describes the most significant direction of data change, and the second principal component describes the less significant direction of data change; after the eigenvalue decomposition of the polarization covariance matrix C2, a set of basis vectors can be obtained according to the eigenvectors, and the basis vectors are projected onto the first principal component and the second principal component, respectively, to obtain the relevant projections of C2 in the two directions; since the principal components constitute an orthogonal basis, the projections of C2 in other principal component directions are regarded as irrelevant projections; the relevant and irrelevant projections are added together to obtain the decomposition result of the polarization covariance matrix C2, that is:

[0065]

[0066] In the formula, C {2，q} represents the value of C2 at position q, C {τ，q} represents the projection of the τth eigenvector of C2 at position q, λ1 and λ2 represent the first and second eigenvalues of C2, respectively, e {1，q} , e {2，q} are the values of the first and second eigenvectors of the matrix C2 at position q, respectively, are the complex conjugate transposes of the first and second eigenvectors of the matrix C2 at position q, respectively, δ1, δ2, ω1, ω2, are constant terms, and I is a 2x2 unit matrix.

[0067] Moreover, the calculation methods of the polarization scattering entropy H, the average scattering angle α and the polarization scattering anti-entropy A in the step 2.3 are as follows:

[0068]

[0069]

[0070]

[0071]

[0072] In the formula, λ1 and λ2 are the first and second eigenvalues of the polarization covariance matrix, respectively, and |e τ (1) | represents the modulus of the first principal polarization component on the τth eigenvector of the polarization covariance matrix.

[0073] Moreover, the cloud region image reconstruction module in step 3.1 includes a local fine texture feature coding submodule and a global attention & receptive field joint enhancement repair submodule. The local fine texture feature coding submodule is a self-encoding model, the encoder is a screening & adjusting down-sampling convolutional layer, the intermediate layer is a residual block with dilated convolution, and the decoder is a screening & adjusting up-sampling convolutional layer. The global attention & receptive field joint enhancement repair submodule is a self-encoding model, including a traditional convolution down-sampling layer, a traditional convolution up-sampling layer, and a high receptive field convolution layer.

[0074] The encoder and the decoder in the local fine texture feature coding submodule selectively transfer useful features through screening & adjusting convolution. The screening & adjusting convolution is performed by a convolution operation unit with a screening mechanism. And a feature extraction unit The two parts are defined as follows:

[0075]

[0076]

[0077]

[0078] In the formula, Φ is an activation function, is an element-level exclusive or operator, are screening & adjusting weights and feature weights in the screening & adjusting convolution, respectively, (y, x) is the coordinate of a specific position, n and m are the size of the screening & adjusting convolution, Ψ is a screening & adjusting function, represents an input image, i.e., the polarized scattering entropy H, the average scattering angle a, and the polarized scattering anti-entropy A extracted in step 2.

[0079] Moreover, in step 3.3, four β k , k∈{0,1,2,3} are initialized to zero, and then the four related feature maps S k , k∈{0,1,2,3} from the local fine texture feature coding submodule are initialized to zero, then they are weighted and summed, and the summation result is input as a new feature map into the next layer network; for the convolutional neural network encoder layer of the global attention & receptive field joint enhancement repair submodule based on Conv-BatchNorm-ReLU, the features are added with an attention mechanism in the original residual term as follows:

[0080] a) The original features are input into the feature extraction function of the convolutional neural network encoder layer to obtain intermediate features

[0081]

[0082] b) Calculate attention coefficients based on intermediate features

[0083]

[0084] c) Weight the intermediate features with attention coefficients to generate new features

[0085]

[0086] d) Perform convolution operation in residual term to get new features

[0087]

[0088] e) Perform batch normalization operation in residual term to get new features

[0089]

[0090] f) Perform ReLU activation operation in residual term to get new features

[0091]

[0092] g) Finally add the residual term to the original feature to get new new features

[0093]

[0094] where, represents a convolution operation on the intermediate feature , and then performs Softmax normalization to get attention coefficients represents the feature obtained by weighting the original feature with attention coefficients.

[0095] Moreover, the high receptive field convolution layer in the global attention & receptive field joint enhancement repair submodule in step 3.4 takes the following formula to represent the consideration of the global structure:

[0096] a) Apply real two-dimensional fast Fourier transform to the input tensor:

[0097]

[0098] Then concatenate the real and imaginary parts:

[0099]

[0100] b) applying a convolution block in the frequency domain:

[0101]

[0102] c) applying an inverse transform to recover the texture structure:

[0103]

[0104]

[0105] constructing a loss function training and optimizing the cloud area image reconstruction module, expressed as:

[0106]

[0107] In the formula, is a global optimization loss function, is a high receptive field perception loss function, and respectively represent the weights of .

[0108] Global optimization loss function The specific formula is as follows:

[0109]

[0110] In the formula, R cloud represents the input cloud detection result, and ⊙ represents element-wise multiplication, respectively represent the real cloud-free image and the predicted cloud-free image;

[0111] High receptive field perception loss function The specific formula is as follows:

[0112]

[0113] In the formula, is an expected value, and ψ hrf represents a pre-trained ResNet50 network with dilated convolution.

[0114] Compared with the prior art, the present application has the following advantages:

[0115] 1) The present application first applies the polarization information of SAR images to the reconstruction of cloud-shading areas of optical remote sensing images, effectively overcoming the problem of rough ground feature details in the repair results obtained by using only the backscattering intensity information of SAR images in the past.

[0116] 2) For the current high-resolution remote sensing image, when image reconstruction is performed using a popular convolutional architecture, the global spatial blur problem caused by the lack of a large enough effective receptive field is solved by adopting a joint attention enhancement and high receptive field perception loss strategy to achieve fine texture reconstruction of the image. BRIEF DESCRIPTION OF DRAWINGS

[0117] Figure 1 Flowchart of the present embodiment.

[0118] Figure 2 Flowchart of the cloud detection module designed by the present application.

[0119] Figure 3 Flowchart of the cloud area image reconstruction module designed by the present application. DETAILED DESCRIPTION

[0120] The present application provides a remote sensing image cloud area reconstruction method based on SAR prior knowledge guidance, which can fully utilize the prior information of SAR image when reconstructing the cloud-occluded area of optical remote sensing image. By using a local fine texture feature coding sub-module and a global attention & receptive field joint enhancement repair sub-module, the consistency repair of the cloud-occluded area in terms of space, texture, details, color, etc. can be realized, and finally a high-fidelity cloud-free image is generated.

[0121] The technical solutions of the present application will be further described below in combination with the drawings and embodiments.

[0122] As shown in Figure 1 , the present application provides a remote sensing image cloud area reconstruction method based on SAR prior knowledge guidance, comprising the following steps:

[0123] Step 1, constructing a cloud detection module, inputting the acquired optical remote sensing image into the cloud detection module for cloud spectrum test and cloud probability calculation to obtain the cloud detection result.

[0124] The optical remote sensing image input into the cloud detection module is first subjected to cloud spectrum test to obtain a preliminary cloud detection result, and then subjected to cloud probability calculation to obtain a cloud edge zone detection result, and then the two are combined to obtain the final cloud detection result, which comprises the following sub-steps:

[0125] Step 1.1, performing cloud spectrum test on the acquired optical remote sensing image to obtain a preliminary cloud detection result.

[0126] The preliminary cloud detection result is R a , r a is a binary image, where 1 represents a cloud pixel and 0 represents a non-cloud pixel. The cloud spectrum test specifically includes basic test, whiteness test, haze detection, rock and desert detection, and cirrus test.

[0127] Step 1.1.1, basic test is performed on the optical remote sensing image.

[0128] A short-wave infrared 2 (SWIR2) band threshold of 0.03, a normalized difference vegetation index (NDVI) threshold of 0.8, and a normalized difference snow index (NDSI) threshold of 0.8 are set respectively to obtain a preliminary cloud detection result R a1 , that is,

[0129]

[0130]

[0131]

[0132] In the formula, p Red , p Green , and p SWIR1 respectively represent reflectivity of a red light band, a green light band, and a short-wave infrared 1 band.

[0133] Step 1.1.2, whiteness test is performed on the optical remote sensing image.

[0134] The whiteness test is performed by calculating a difference value of each visible light band from an average value of the visible light band divided by the average value of the visible light band, that is,

[0135]

[0136] meanVis = (p Red + p Green + p Blue ) / 3 (5)

[0137] In the formula, whiteness represents the whiteness test, meanVis represents the average value of each visible light band, p Red , p Green , and p Blue are reflectivity of a red band, a green band, and a blue light band respectively.

[0138] A threshold value of the result of the whiteness test is set to 0.7 to obtain a preliminary cloud detection result R a2 = whiteness < 0.7.

[0139] Step 1.1.3, haze detection is performed on the optical remote sensing image.

[0140] A score greater than 0 is set as a cloud pixel by a haze optimized transformation (HOT) algorithm to obtain a preliminary cloud detection result R a3 = HOT > 0, and a specific calculation formula of the HOT is as follows:

[0141] HOT = p Blue -0.5 x p Red -0.08 (6)

[0142] wherein p Red , p Blue are reflectance of red band, blue band respectively.

[0143] Step 1.1.4, rock desert detection is performed on the optical remote sensing image.

[0144] By calculating the ratio of the near-infrared (NIR) band and SWIR1 of the image, the pixels with the result greater than 0.75 are set as cloud pixels, and the preliminary cloud detection result R a4 is obtained, that is:

[0145] R a4 = p NIR / p SWIR1 > 0.75 (7)

[0146] wherein p NIR , p SWIR1 are reflectance of near-infrared band, SWIR1 band respectively.

[0147] Step 1.1.5, cirrus cloud test is performed on the optical remote sensing image.

[0148] The reflectance of the cirrus cloud band is divided by 0.04, and the pixels with the result greater than 0.01 are set as cloud pixels, and the preliminary cloud detection result R a5 is obtained, that is:

[0149] R a5 = p Cirrus / 0.04 > 0.01 (8)

[0150] wherein p Cirrus represents reflectance of cirrus cloud band.

[0151] Step 1.1.6, the five preliminary cloud detection results obtained in steps S1.1.1-S1.1.5 are fused to obtain the final preliminary cloud detection result.

[0152] The final preliminary cloud detection result R a is represented as:

[0153] R a = R a1 ∪ R a2 ∪ R a3 ∪ R a4 ∪ R a5 (9)

[0154] Step 1.1.7, calculate the difference between cloud and land pixels in spectral reflectance to obtain the spectral change probability.

[0155] Spectral change probability L vari is calculated as follows:

[0156] L vari = 1 - max(|NDVI|, |NDSI|, |whitenessD (10)

[0157] where NDVI is the normalized vegetation index, NDSI is the normalized snow index, and whiteness is the whiteness test result.

[0158] Step 1.1.8, normalize the HOT result obtained in step 1.1.3.

[0159] The normalization calculation formula is as follows:

[0160]

[0161] where HOT low and HOT high represent the 17.5% quantile and 82.5% quantile of the HOT value, respectively.

[0162] Step 1.2, calculate the cloud probability based on the preliminary cloud detection result, the spectral change probability, and the HOT normalization result to obtain the cloud edge zone detection result.

[0163] Combine the preliminary cloud detection result R a5 obtained in step 1.1.5, the spectral change probability L vari obtained in step 1.1.7, and the HOT normalization result L HOT obtained in step 1.1.8 to obtain the cloud probability Cloud p Set the value of Cloud p greater than 0.8 as a cloud pixel to obtain the cloud edge zone detection result R b .

[0164] The cloud edge zone detection result R b is a binary image, where 1 represents a cloud pixel and 0 represents a non-cloud pixel. R b is calculated as follows:

[0165] R b = Cloud o = L vari × L HOT + R a5 × 0.5 > 0.8 (12)

[0166] Step 1.3, fuse the preliminary cloud detection result R

[0167] with the cloud edge zone detection result R a obtained in step 1.2, to obtain the final optical remote sensing image cloud detection result R b . cloud R cloud is a binary image, where 1 represents a cloud pixel, and 0 represents a non-cloud pixel. R cloud is as follows:

[0168] R cloud = R a ∩ R b (13)

[0169] Step 2, construct a SAR feature extraction module, input the acquired near-time SAR image into the SAR feature extraction module, construct a polarization covariance matrix of the SAR image, and extract relevant polarization scattering features.

[0170] The extracted polarization scattering features include polarization scattering entropy H, average scattering angle a, and polarization scattering anti-entropy A.

[0171] Step 2.1, construct a polarization covariance matrix of the SAR image.

[0172] The polarization covariance matrix C2 is constructed as follows:

[0173]

[0174] In the formula, V represents the vertical polarization state of the SAR antenna; H represents the horizontal polarization state of the SAR antenna; S VH represents the relationship between the incident field and the scattering field in the polarization electromagnetic wave state of transmitting a vertical signal and receiving a horizontal signal; S VV represents the relationship between the incident field and the scattering field in the polarization electromagnetic wave state of transmitting a vertical signal and receiving a vertical signal; * represents complex conjugate; <·> represents statistical average; |S VH | and |S VV | represent the phase information of the SAR image under different polarization states, and i is an imaginary unit; represents the polarization information of the SAR image.

[0175] Step 2.2, decompose the polarization covariance matrix C2 into the sum of relevant projections and irrelevant projections in different directions.

[0176] The eigenvalue decomposition of the polarimetric covariance matrix C2 gives the eigenvectors and eigenvalues. The eigenvectors are unit vectors that describe the variation of the data in different directions, while the eigenvalues indicate the variation degree of the data in the corresponding eigenvector direction. The eigenvectors can be sorted according to the corresponding eigenvalues, and a set of principal components can be obtained. Among these principal components, the first principal component describes the most significant direction of data variation, the second principal component describes the less significant direction of data variation, and so on. After the eigenvalue decomposition of the polarimetric covariance matrix C2, a set of basis vectors can be obtained according to the eigenvectors. The projection of these basis vectors onto the first and second principal components gives the correlation projections of C2 in these two directions. At the same time, since the principal components constitute an orthogonal basis, the projections of C2 in other principal component directions can be regarded as uncorrelated projections. Adding these correlated and uncorrelated projections together gives the decomposition result of the polarimetric covariance matrix C2, that is:

[0177]

[0178] where C {2，q} represents the value of C2 at position q, C {τ，q} represents the projection of the τth eigenvector of C2 at position q, λ1 and λ2 represent the first and second eigenvalues of C2, respectively, e {1，q} and e {2，q} are the values of the first and second eigenvectors of matrix C2 at position q, respectively, are the complex conjugate transposes of the first and second eigenvectors of matrix C2 at position q, respectively, δ1, δ2, ω1, ω2, are constant terms, and I is a 2x2 identity matrix, that is:

[0179]

[0180] Step 2.3, using the decomposition result of the polarimetric covariance matrix, the polarimetric scattering entropy, the average scattering angle and the polarimetric scattering anti-entropy are calculated.

[0181] The calculation methods of the polarimetric scattering entropy H, the average scattering angle α and the polarimetric scattering anti-entropy A are as follows:

[0182]

[0183]

[0184]

[0185]

[0186] where λ1 and λ2 are the first and second eigenvalues of the polarimetric covariance matrix, respectively, and |eτ (1) | represents the modulus value of the first principal polarized component of the polarization covariance matrix on the τth eigenvector.

[0187] Step 3, constructing a cloud area image reconstruction module, inputting the acquired optical remote sensing image, the cloud detection result obtained in step 1 and the polarization scattering feature obtained in step 2 into the cloud area image reconstruction module in turn, realizing spatial structure feature reconstruction and global consistency feature repair of the optical remote sensing image cloud shadow area, and finally generating a high-fidelity cloud-free image.

[0188] The cloud area image reconstruction module includes a local fine texture feature coding submodule and a global attention & receptive field joint enhancement repair submodule. The local fine texture feature coding submodule is a self-encoding model, including 3 screening & adjusting down-sampling convolutional layers (encoder), 3 layers of residual blocks with dilated convolution (intermediate layer) and 3 screening & adjusting up-sampling convolutional layers (decoder). The global attention & receptive field joint enhancement repair submodule is a self-encoding model, having 3 traditional convolutional down-sampling layers, 3 traditional convolutional up-sampling layers, and 9 high-receptive field convolutional layers.

[0189] Step 3.1, using the local fine texture feature coding submodule to map the polarization scattering entropy, average scattering angle and polarization scattering anti-entropy extracted in step 2 to the feature space.

[0190] For the encoder and decoder in the local fine texture feature coding submodule, the screening & adjusting convolution selectively transfers useful features. The screening & adjusting convolution is composed of a convolution operation unit with a screening mechanism and a feature extraction unit , which are defined as follows:

[0191]

[0192]

[0193]

[0194] In the formula, Φ is an activation function, is an element-level exclusive or operator, are the screening & adjusting weights and feature weights in the screening & adjusting convolution respectively, (y, x) is the coordinate of a specific position, n and m are the size of the screening & adjusting convolution, Ψ is a screening & adjusting function, and the ReLU function is selected in this embodiment, represents an input image, and corresponds to the polarization scattering entropy H, the average scattering angle α and the polarization scattering anti-entropy A extracted in step 2.

[0195] Step 3.2, the last 4 feature maps from coarse to fine S k , k e {0, 1, 2, 3} are sequentially output from the last intermediate layer of the local fine texture feature encoding submodule and 3 decoder layers.

[0196] Step 3.3, based on zero initialization residual addition feature fusion technology, gradually add the 4 feature maps from coarse to fine S k , k e {0, 1, 2, 3} obtained in step 3.2 to the global attention & receptive field joint enhancement repair submodule.

[0197] Specifically, 4 zero initialized β k , k e {0, 1, 2, 3} are used to initialize the 4 related feature maps S k , k e {0, 1, 2, 3} from the local fine texture feature encoding submodule with zero tensor, then they are weighted sum, and the sum result is input as a new feature map to the next layer network. For the features of the convolutional neural network encoder layer in the global attention & receptive field joint enhancement repair submodule based on Conv-BatchNorm-ReLU , the attention mechanism is added to the original residual term as follows:

[0198] a) The original feature is subjected to the feature extraction function of the convolutional neural network encoder layer to obtain the intermediate feature

[0199]

[0200] b) Calculate the attention coefficient based on the intermediate feature

[0201]

[0202] c) Use the attention coefficient to weight the intermediate feature to generate a new feature

[0203]

[0204] d) Convolution operation in the residual term to get new features

[0205]

[0206] e) Batch normalization operation in the residual term to get new features

[0207]

[0208] f) ReLU activation operation is performed on the residual term to obtain new features

[0209]

[0210] g) Finally, the residual term is added to the original features to obtain new new features

[0211]

[0212] wherein, represents a convolution operation on the intermediate features , and then a Softmax normalization is performed to obtain attention coefficients represents the features obtained by weighting the original features using the attention coefficients; using The network can enhance the learning and use of important features, and increase the expression ability and complexity of the network.

[0213] Step 3.4, input the cloud-containing optical remote sensing image and the cloud detection result obtained in step 1 into the global attention & receptive field joint enhancement repair sub-module to realize repair of the cloud-shielded area and other areas not shielded by clouds in terms of global consistency, and obtain the final cloud-removed image.

[0214] The repair of the cloud-shielded area and other areas not shielded by clouds in terms of global consistency includes spatial structure features, color consistency features, etc. The global attention & receptive field joint enhancement repair sub-module is an auto-encoding network model, which has 3 traditional convolution down-sampling layers, 3 traditional convolution up-sampling layers, and 9 high receptive field convolution layers. The construction of the high receptive field convolution layer mainly utilizes the fast Fourier transform, which has a receptive field covering the height (H), width (W), and channel number (C) of the image, which means that this layer can effectively capture the relationship between different positions and channels in the image, thereby extracting more representative features. Specifically, the high receptive field convolution layer adopts the following formula to represent the consideration of global structure:

[0215] a) Apply real two-dimensional fast Fourier transform (Real FFT2d) to the input tensor:

[0216]

[0217] Then, the real and imaginary parts are spliced (ComplesToReal):

[0218]

[0219] b) applying a convolution block (ReLU°BN°Conv1x1) in the frequency domain:

[0220]

[0221] c) applying an inverse transform to recover the texture structure (RealToComples, Inverse Real FFT2d):

[0222]

[0223]

[0224] Constructing the loss function The cloud area image reconstruction module is trained and optimized, expressed as:

[0225]

[0226] In the formula, R is a global optimization loss function, is a high receptive field perception loss function, and respectively represent the weights.

[0227] Global optimization loss function The specific formula is as follows:

[0228]

[0229] In the formula, R cloud represents the input cloud detection result, and represents element-by-element multiplication, respectively represent the real cloud-free image and the predicted cloud-free image.

[0230] High receptive field perception loss function The specific formula is as follows:

[0231]

[0232] In the formula, is an expected value, and ψ hrf represents a pre-trained ResNet50 network with dilated convolution.

[0233] The specific embodiments described herein are merely illustrative of the spirit of the present application. Those skilled in the art of the present application can make various modifications or supplements to the described specific embodiments or use similar ways to replace them, without deviating from the spirit of the present application or exceeding the scope defined by the appended claims.

Claims

1. A method for cloud region reconstruction of remote sensing image based on SAR prior knowledge guidance, characterized in that, The method comprises the following steps: Step 1, constructing a cloud detection module, inputting the obtained optical remote sensing image into the cloud detection module for cloud spectrum testing and cloud probability calculation to obtain a cloud detection result; Step 1.1, performing cloud spectrum testing on the obtained optical remote sensing image to obtain a preliminary cloud detection result; The cloud spectrum testing comprises basic testing, whiteness testing, haze detection, rock and desert detection and cirrus cloud testing; Step 1.2, extracting a cloud edge zone in the image by calculating a cloud probability to obtain a cloud edge zone detection result; Step 1.3, fusing the preliminary cloud detection result and the cloud edge zone detection result to obtain a final optical remote sensing image cloud detection result; Step 2, constructing a SAR feature extraction module, inputting the obtained near-time SAR image into the SAR feature extraction module, constructing a polarized covariance matrix of the SAR image and extracting relevant polarized scattering features; The extracted polarimetric scattering features include polarimetric scattering entropy , mean scattering angle , and polarimetric scattering negentropy , including the following steps: Step 2.1, constructing a polarized covariance matrix of the SAR image; Polarization covariance matrix The construction is as follows: （14） where V represents the vertical polarization state of the SAR antenna; H represents the horizontal polarization state of the SAR antenna; represents the relationship between the incident field and the scattering field in the state of transmitting a vertical signal and receiving a horizontal signal polarized electromagnetic wave; represents the relationship between the incident field and the scattering field in the state of transmitting a vertical signal and receiving a vertical signal polarized electromagnetic wave; represents a complex conjugate; represents a statistical average; , represents the phase information of the SAR image in different polarization states, and ; is an imaginary unit; represents the polarization information of the SAR image; Step 2.2, decomposing the polarized covariance matrix into a sum of relevant projections and irrelevant projections in different directions; Polarization covariance matrix Eigenvalue decomposition is performed to obtain eigenvectors and eigenvalues, and the eigenvectors are sorted according to the corresponding eigenvalues to obtain a set of principal components, wherein the first principal component describes the most significant direction of data change, and the second principal component describes the less significant direction of data change; the polarization covariance matrix After eigenvalue decomposition, a set of basis vectors can be obtained according to the eigenvectors, and the basis vectors are projected onto the directions of the first principal component and the second principal component, respectively, to obtain Correlated projections in the two directions; since the principal components constitute an orthogonal basis, the projections of the basis vectors in the directions of other principal components are regarded as uncorrelated projections; the correlated and uncorrelated projections are added together to obtain the decomposition result of the polarization covariance matrix , that is: （15） wherein denotes the value of the at position denotes the projection of the Step 2.3, calculating polarized scattering entropy, average scattering angle and polarized scattering anti-entropy by using the polarized covariance matrix decomposition result; th eigenvector of the matrix at position denotes the first, second eigenvector of the matrix , respectively, are the values of the first, second eigenvector of the matrix at position are the complex conjugate transposes of the first, second eigenvector of the matrix at position , , , , , , is a constant term, is a 2x2 identity matrix; Step 3, constructing a cloud area image reconstruction module, inputting the obtained optical remote sensing image, the cloud detection result obtained in step 1 and the polarized scattering features obtained in step 2 into the cloud area image reconstruction module in sequence to realize spatial structure feature reconstruction and global consistency feature repair of the optical remote sensing image cloud shadow area and finally generate a high-fidelity cloud-free image; The cloud area image reconstruction module comprises a local fine texture feature coding submodule and a global attention & receptive field joint enhancement repair submodule; the local fine texture feature coding submodule is a self-encoding model, the encoder is a screening & adjusting down-sampling convolutional layer, the intermediate layer is a residual block with dilated convolution and the decoder is a screening & adjusting up-sampling convolutional layer; the global attention & receptive field joint enhancement repair submodule is a self-encoding model, comprising a traditional convolution down-sampling layer, a traditional convolution up-sampling layer and a high receptive field convolution layer. Step 1.1 specifically comprises the following sub-steps: 2.The SAR prior knowledge guided remote sensing image cloud area reconstruction method according to claim 1, characterized in that: Step 1.1.1, performing basic testing on the optical remote sensing image; Step 1.1.2, performing whiteness testing on the optical remote sensing image; Let the initial cloud detection results be , is a binary image, where 1 represents cloud pixels, and 0 represents non-cloud pixels, respectively set the shortwave infrared 2 band threshold value to , the normalized vegetation index threshold value to , and the normalized snow cover index threshold value to , to obtain the initial cloud detection results of the optical remote sensing image , that is: （1） （2） （3） In the formula, , , respectively represent the reflectivity of red light band, green light band, short wave infrared 2 band, NDVI is normalized vegetation index, and NDSI is normalized snow cover index. The whiteness testing is performed by calculating a difference value of each visible light band and the average value of the visible light band divided by the average value of each visible light band, that is: Step 1.1.3, performing haze detection on the optical remote sensing image; （4） （5） wherein whiteness represents the whiteness test, meanVis represents the average of the various visible light bands, respectively the reflectance of the red, green, and blue light bands. The result of the whiteness test is thresholded to obtain a preliminary cloud detection result ; Step 1.1.4, performing rock and desert detection on the optical remote sensing image; The score greater than 0 is set as a cloud pixel by the haze optimization conversion HOT algorithm to obtain a preliminary cloud detection result The specific calculation formula of HOT is as follows: （6） In the formula, respectively, the reflectivity of red band, blue light band; Step 1.1.5, performing cirrus cloud testing on the optical remote sensing image; The cloud pixels are obtained by calculating the ratio of the near-infrared band and the short-wave infrared 1 band of the image, and the obtained result greater than The preliminary cloud detection result is obtained by setting the pixel as a cloud pixel That is: （7） In the formula, , are reflectance of near-infrared band and short-wave infrared 1 band, respectively. Step 1.1.6, fusing the five preliminary cloud detection results obtained in steps 1.1.1-1.1.5 to obtain a final preliminary cloud detection result; Divide the reflectivity of the cirrus band by 0.04, and obtain a result greater than The pixel is set as a cloud pixel, and a preliminary cloud detection result is obtained That is （8） wherein Rcindicates the reflectivity of the cirrus cloud band; Step 1.1.7, calculating a difference between cloud and land pixels in spectral reflectance to obtain a spectral change probability; Final preliminary cloud detection results is represented as: （9） In the formula, NDVI is a normalized vegetation index, NDSI is a normalized snow index and whiteness is a whiteness test result; Spectral change probability The calculation formula is as follows: （10） ​ Step 1.1.8, normalize the HOT result obtained in step 1.1.3; The normalization calculation formula is as follows: （11） wherein represents the normalized HOT result, and represents the HOT value 17.5% quantile and 82.5% quantile, respectively. 3.The SAR prior knowledge guided remote sensing image cloud area reconstruction method according to claim 2, characterized in that: Step 1.2 is to combine the preliminary cloud detection result obtained in step 1.1.5, the spectral change probability index obtained in step 1.1.7 and the HOT normalization result obtained in step 1.1.8 to obtain a cloud probability , the value greater than is set as a cloud pixel, and a cloud edge zone detection result is obtained ; and ; Cloud edge zone detection results is a binary image where 1 represents a cloud pixel and 0 represents a non-cloud pixel, The calculation is as follows: （12） wherein is a preliminary cloud detection result, is a probability of spectral change, is a HOT normalization result.

4. The SAR prior knowledge guided remote sensing image cloud region reconstruction method of claim 3, wherein: Step 1.3 is to fuse the preliminary cloud detection result obtained in step 1.1.5 with the cloud edge zone detection result obtained in step 1.2 to obtain a final cloud detection result of the optical remote sensing image , is a binary image, wherein 1 represents a cloud pixel, and 0 represents a non-cloud pixel, The calculation is as follows: (13) wherein, is a preliminary cloud detection result, is a cloud edge zone detection result.

5. The SAR prior knowledge guided remote sensing image cloud region reconstruction method of claim 1, wherein: The polarized scattering entropy in step 2.3 The average scattering angle and the polarized scattering anti-entropy are calculated as follows: （17） （18） （19） （20） wherein are the first and second eigenvalues of the polarization covariance matrix, respectively, denotes the modulus of the first principal polarized component of the polarization covariance matrix on the τ first eigenvector.

6. The SAR prior knowledge guided remote sensing image cloud region reconstruction method of claim 5, wherein: Step 3 includes the following steps: Step 3.1, using the local fine texture feature coding submodule to map the polarization scattering entropy, average scattering angle and polarization scattering anti-entropy extracted in step 2 into the feature space; The encoders and decoders in the local fine texture feature encoding sub-module construct a screening & conditioning convolution that selectively transfers useful features, the screening & conditioning convolution from a convolution operation unit with a screening mechanism and a feature extraction unit consists of two parts, which are defined as follows: （21） （22） （23） wherein, is an activation function, is an element-wise XOR operator, , are the screening & conditioning weights and feature weights in the screening & conditioning convolution, respectively, is the coordinate of a specific position, , is the size of the screening & conditioning convolution, is a screening & conditioning function, represents the input image, i.e. the polarized scattering entropy extracted in step 2 , the average scattering angle and the polarized scattering de-entropy ; Step 3.2, the last 4 feature maps from coarse to fine are output from the last intermediate layer and 3 decoder layers of the local fine texture feature encoding submodule in turn k e {0, 1, 2, 3} Step 3.3, based on zero initialization residual addition feature fusion technology, gradually add the four feature maps from coarse to fine obtained in step 3.2 k∈{0,1,2,3} to the global attention & receptive field joint enhancement repair submodule; Step 3.4, input the cloud-containing optical remote sensing image and the cloud detection result obtained in step 1 into the global attention & receptive field joint enhancement repair submodule to realize the repair of the cloud-shielded area and other areas not shielded by the cloud in the global consistency, and obtain the final cloud-removed image.

7. The SAR prior knowledge guided remote sensing image cloud region reconstruction method of claim 6, wherein: The four zero-initialized , k∈{0,1,2,3} four relevant feature maps from the local fine texture feature encoding sub-module k∈{0,1,2,3} are initialized with all-zero tensors, then they are weighted summed, and the summation result is input into the next layer network as a new feature map; for the convolutional neural network encoder layer in the global attention & receptive field joint enhancement repair sub-module based on Conv-BatchNorm-ReLU The features of the original residual term are added with attention mechanism in the following way: a) original features feature extraction function through a convolutional neural network encoder layer obtaining intermediate features : （24） b) calculating attention coefficients based on the intermediate features : （25） c) weighting the intermediate features with attention coefficients to generate new features : （26） d) performing a convolution operation in the residual term to obtain new features : （27） e) performing a batch normalization operation in the residual term to obtain new features : （28） f) performing a ReLU activation operation in the residual term to obtain new features : （29） g) finally adding the residual term to the original feature to obtain a new new feature : （30） In the formula, is based on the intermediate feature A convolution operation is performed once, and then a Softmax normalization is performed to obtain the attention coefficient ; The feature obtained by weighting the original feature with the attention coefficient is represented.

8. The SAR prior knowledge guided remote sensing image cloud region reconstruction method of claim 7, wherein: The high receptive field convolution layer in the global attention & receptive field joint enhancement repair submodule in step 3.4 takes the following formula to express the consideration of the global structure: a) Apply real two-dimensional fast Fourier transform to the input tensor: （31） Then the real part and the imaginary part are spliced together: （32） b) Apply convolution block in frequency domain: （33） c) Apply inverse transform to restore texture structure: （34） （35） Constructing loss function Training and optimizing the cloud area image reconstruction module, expressed as: （36） wherein is a global optimization loss function, is a high-receptive field perception loss function, and respectively represent , the weight of Global optimization loss function The specific formula is as follows: （37） In the formula, represents the input cloud detection result, represents element-wise multiplication, , respectively represent the real cloud-free image and the predicted cloud-free image; High receptive field perception loss function The specific formula is as follows: （38） In the formula, is the expected value, represents a pre-trained ResNet50 network with dilated convolution.

Images