CN103049898B - Method for fusing multispectral and full-color images with light cloud - Google Patents

Abstract

The invention provides a method for fusing multispectral and full-color images with light cloud. The problem that the cloud and mist area is interfered by a cloud layer after the multispectral and full-color images with light cloud are fused is mainly solved. The method comprises the following implementation steps of: sampling and filtering the multispectral image with light cloud, and obtaining background images of the multispectral and full-color images; respectively removing light cloud of the multispectral and full-color images with light cloud; performing PCA conversion on the multispectral image in which the cloud is removed, and performing Shearlet decomposition on the converted first principal component image and the full-color image; taking a low-frequency coefficient of the first principal component as a low-frequency coefficient of a fusion component, and taking a high-frequency coefficient of the full-color image as a high-frequency coefficient of the fusion component; and performing reversible PCA conversion on the fusion component and the other components subjected to the PCA conversion, and obtaining the fusion image. The method has the advantages of high definition of light cloud area of the fusion image and high spectrum retainability and can be used for military target recognition, weather and environment monitoring, land utilization, urban planning and prevention and reduction of natural disasters.

Description

Multispectral and Panchromatic Image Fusion Method with Thin Clouds

technical field

The invention belongs to the field of intelligent image processing, relates to thin cloud removal and image fusion methods, and can be used in technologies in multiple fields such as military target recognition, meteorological monitoring, environmental monitoring, land utilization, urban planning, and disaster prevention and mitigation.

Background technique

Since the United States implemented the Earth Resources Satellite Program in 1972, satellite technology has developed rapidly around the world. So far, major developed countries and a few developing countries in the world, including China and India, have successively launched hundreds of satellites. The operating band covers Visible light to invisible near-infrared, short-wave infrared, mid-infrared, far-infrared, microwave and other broad frequency domains. These satellites transmit massive amounts of satellite remote sensing data covering the globe to satellite ground stations and mobile receiving stations scattered around the world every day. Due to the spatial continuity and time sequence of satellite remote sensing data, so far, it is the only means that can provide dynamic observation data on a global scale.

At present, most of the widely used satellite remote sensing data are optical images, although optical images generally have the characteristics of large amount of information, high resolution and stable images. But at the same time, due to the common presence of clouds and fog in the atmosphere during the imaging process, remote sensing satellites often record the image data together with cloud and fog information. Due to the occlusion of clouds and fog, the sensor cannot obtain the information of ground objects in the area covered by clouds and fog, which seriously affects the quality of remote sensing optical images. When part of the image is covered by thick clouds and fog, the information of ground objects cannot be received by the sensor. Relatively thin clouds, the sensor can still receive part of the ground object information, but this incomplete information has a serious impact on the recognition and classification of the target in the image, and the accuracy of the ground object information extraction.

In order to effectively improve the utilization of remote sensing optical images, there are many methods on the market to reduce or remove the influence of thin clouds and haze on images. At present, the main methods used in the market are: multi-spectral image method, multi-image overlay algorithm, method based on image enhancement, and method based on image restoration.

The multi-spectral method uses a special sensor, or some bands in the multi-spectral image have strong sensitivity to clouds and fog to detect cloud and fog information. Then the cloud and fog information is subtracted from the original image to obtain the cloud-removed image. This method can effectively remove the cloud and fog in the image, but it needs to add sensors or wave bands that are sensitive to the cloud and fog, and the cost is high, which limits the application of this method.

The multi-image superposition method is to obtain information images by superimposing images taken in different seasons and at different times in the same area. Since the ground object information in the same area often changes at different times, it seriously affects the interpretation of the superimposed image.

The cloud and fog removal method based on image restoration is to find out the corresponding inverse process by analyzing the mechanism and process of thin cloud degradation image, so as to obtain the original image. However, to process images from the perspective of image restoration, one must be familiar with the mechanism and process of thin cloud degraded images. Since the degradation degree of thin cloud degraded images is related to the distance between the target and the camera, relevant auxiliary information, such as the atmospheric extinction coefficient, must be combined during processing. , target and camera distance, etc., and the acquisition cost of these information is high, and it is difficult to be widely used.

There are usually two methods for removing clouds and fog based on image enhancement, homomorphic filtering and low-frequency filtering. Homomorphic filtering is a method that compresses the dynamic range of the image and enhances the high-frequency components of the image to remove thin clouds. This method does not distinguish the details of the image, and uses a single filter for enhancement, which has poor adaptability; the low-frequency filtering method performs Gaussian low-pass filtering on the image to obtain the background image of the image, and subtracts the original image from the background image to obtain To achieve the purpose of removing thin clouds. This method has a certain removal effect on thin clouds and mist, but it also weakens the background of the image. Although the compensation method was used to improve it, from the results, this method still produces information on the image after removing clouds and fog. loss.

Since the image acquired by a single sensor is difficult to meet the actual needs in terms of spectral information and resolution, it is necessary to use the redundancy of sensor image data to fuse the image information of different sensors into one image through fusion technology. , such as multispectral and panchromatic images, multispectral images have better spectral information, but their spatial resolution is low; panchromatic images have higher spatial resolution, but their spectral information is not rich enough, so fusion technology is needed , to obtain high-resolution multispectral images to obtain a more accurate description of surface object information. However, when there is cloud and fog information in multispectral and panchromatic images, the fusion technology cannot remove the cloud and fog information in the image. For the fusion of multi-spectral and panchromatic images with thin clouds, the current market is to firstly remove the clouds and fog of these two images, and then perform fusion processing. When removing clouds and fog from the images acquired by a single sensor, from the perspective of the cost of acquiring image information, the method of image enhancement should be selected, but the loss of image information caused by image enhancement is inevitable. Therefore, using the redundancy between multispectral and panchromatic images to find a more effective fusion method for thin cloud images is an urgent problem in the market.

Contents of the invention

The purpose of the present invention is to address the deficiencies in the fusion of images with thin clouds in the above-mentioned prior art, and propose a multi-spectral and panchromatic image fusion method with thin clouds, so as to reduce the impact of clouds on thin clouds after the fusion of multi-spectral and panchromatic images. The interference of regional information reduces the loss of fused image information, improves the clarity of fused images, and more accurately describes ground object information.

To achieve the above object, the main contents of the present invention are as follows:

1) Upsampling the multispectral image with thin clouds, the size of the sampled multispectral image is the same as the panchromatic image with thin clouds;

2) Down-sampling, Gaussian low-pass filtering and up-sampling are respectively performed on the sampled multispectral and panchromatic images to obtain the background image B ₁ of the multispectral image and the background image B ₂ of the panchromatic image, in which The number of bands of the background image is the same as that of the multispectral image;

3) The thin cloud removal operation is performed on each band of the multispectral image and the panchromatic image, and the multispectral image I ₁ with thin cloud removed and the panchromatic image I ₂ with thin cloud removed are obtained;

4) Perform PCA transformation on the multispectral image I ₁ after removing thin clouds to obtain each component image in Represents the i-th principal component of the multispectral image after PCA transformation, i=1, 2,..., n, n is the total number of components;

5) For the first principal component image obtained after PCA transformation and the panchromatic image I ₂ with thin clouds removed, respectively decompose by Shearlet transform, and the first principal component image Decompose into a low-frequency coefficient x ₁ and multiple directional subband coefficients y ₁ , y ₂ , ..., y _m , decompose the panchromatic image I ₂ with thin clouds removed into a low-frequency coefficient x ₂ and multiple directional sub-bands Coefficients z ₁ , z ₂ , ... z _m ;

6) For the background image of each band of the multispectral image obtained in step 2), calculate the average gray value to obtain the composite background image of the multispectral image

7) Synthetic background image based on multispectral image Establish a weight matrix w ₁ , and establish a weight matrix w ₂ according to the background image B ₂ of the panchromatic image;

8) For each directional subband coefficient z ₁ , z ₂ , ..., z _m of the panchromatic image I ₂ , multiply the weight matrix w ₁ and w ₂ to obtain the directional subband of the fused first component image Coefficients l ₂ , l ₂ ,..., l _m , using the low-frequency coefficient x ₁ of the first principal component image as the low-frequency coefficient k of the first principal component image after fusion;

9) Perform inverse Shearlet transform on the low-frequency coefficient k and multiple directional subband coefficients l ₁ , l ₂ , ..., l _m of the fused first principal component image to obtain the fused first principal component image

10) The fused first principal component image and other component images obtained in step 4) except the first principal component A new data set is formed, and the inverse PCA transformation is performed on the new data set to obtain the fused image I.

Compared with the prior art, the present invention has the following advantages:

(a) Since the present invention uses the up-sampled image to obtain the background image, the approximate background image of the original image is obtained by sampling the up-sampled image, and the information of the up-sampled image is very smooth, and the approximate background image obtained by sampling The image can better approximate the characteristics of the real background image, overcome the problem of high computational time complexity in the filtering process of the traditional algorithm, and thus improve the speed of removing thin clouds from the image.

(b) Since the present invention uses the background image to establish a weight matrix to enhance the detailed information of the thin cloud area of the fusion image, and utilizes the characteristics of the thin cloud area in the image to have a higher gray value, it overcomes the thin cloud area in the traditional fusion image. The loss of information leads to the blurring of the fused image, thereby improving the clarity of the fused image.

Description of drawings

Fig. 1 is a flow chart of the multispectral and panchromatic image fusion method with thin clouds in the present invention;

Fig. 2 is the QuickBird satellite image with thin cloud of the present invention;

Fig. 3 is the fused image obtained by the method of the present invention.

detailed description

With reference to Fig. 1, the multispectral and panchromatic image fusion method steps that have thin cloud of the present invention are as follows:

Step 1. Upsampling the multispectral image with thin clouds. The size of the upsampled multispectral image is the same as that of the panchromatic image with thin clouds. The original multispectral image is shown in Figure 2(a) As shown, denoted as Y ₁ , the original full-color image is shown in Figure 2(b), denoted as Y ₂ .

Step 2: Sampling and Gaussian low-pass filtering are performed on the upsampled multispectral image and the original panchromatic image respectively to obtain the background image B ₁ of the multispectral image with thin clouds and the background image B ₂ of the panchromatic image.

2.1) Downsampling is performed on the upsampled multispectral image and the original panchromatic image respectively, and the size of the downsampled multispectral and panchromatic image is 1/3 of the size of the original panchromatic image;

2.2) Gaussian low-pass filtering is performed on the downsampled multispectral image and panchromatic image, respectively, to obtain the background image of the downsampled multispectral image and panchromatic image, and the background image of the downsampled multispectral image The number of bands is the same as that of the multispectral image;

2.3) Up-sampling the down-sampled background image to obtain the background image B ₁ of the multispectral image and the background image B ₂ of the panchromatic image, and the size of the obtained images is the same as that of the original panchromatic image.

Step 3, each band of the original multispectral image _Y1 and the original panchromatic image are respectively subjected to thin cloud removal operations.

3.1) Calculate the gray mean value of each band image of the original multispectral image Y ₁ and the gray mean value of the original panchromatic image Y ₂ , that is, sum the gray value of all pixels in the image, and then divide it by the total number of pixels in the image ;

3.2) For each band of the original multispectral image _Y1 , subtract the corresponding band of the background image B1, and add the gray mean value _of the corresponding band to obtain the multispectral image _I1 without thin clouds;

3.3) For the original panchromatic image Y ₂ , subtract its background image B ₂ and add the mean value of its gray level to obtain a panchromatic image I ₂ without thin clouds.

Step 4, perform PCA transformation on the multispectral image I ₁ after thin clouds are removed, and obtain each component image in Indicates the ith principal component of the multispectral image after PCA transformation, i=1, 2,..., n, n is the total number of components.

Step 5, the first principal component image obtained after PCA transformation and the panchromatic image I ₂ with thin clouds removed, respectively decompose by Shearlet transform, and the first principal component image Decompose into a low-frequency coefficient x ₁ and m direction sub-band coefficients y ₁ , y ₂ , ..., y _m , decompose the panchromatic image I ₂ with thin clouds removed into a low-frequency coefficient x ₂ and m direction sub-bands Coefficients z ₁ , z ₂ , . . . , z _m .

Shearlet transform is one of the new generation of multi-scale geometric analysis tools. The specific decomposition method is as follows:

5.1) Scale-decompose the image through the non-subsampling Laplacian pyramid transform, and respectively divide the first principal component image And the panchromatic image _I2 is decomposed into a low-frequency coefficient and multiple high-frequency coefficients, the number of high-frequency coefficients is the same as the number of decomposition layers, and the number of layers of scale decomposition is not limited. In order to make the decomposed low-frequency coefficients and panchromatic The image has a higher correlation. In this example, the number of decomposition layers is determined to be 4;

5.2) Decompose each high-frequency coefficient into multiple directions through the Shear filter, and the choice of the number of directions is not limited. Considering the time complexity of the calculation and the needs of the algorithm, this example determines from the coarse scale to the fine scale. The 4 high-frequency coefficients are decomposed into 6, 6, 10 and 10 directions respectively, and a total of 22 direction subband coefficients are obtained.

Step 6, according to the background image B ₁ of the multispectral image in step 2, calculate the gray mean value of all bands to obtain the composite background image of the multispectral image

The synthetic background image is obtained by the following formula:

${\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{<span\; class="notranslate">\; \&prime;</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; q</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\frac{\underset{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; N</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; q</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; N</span>}}$

in Indicates the gray value of the k- _th band of the background image B1 at coordinates (p, q), Represents the composited background image Gray value at coordinates (p,q), N is the number of bands of the multispectral image.

Step 7, based on the synthetic background image of the obtained multispectral image Create a composite background image The weight matrix w ₁ of .

7.1) Create a synthetic background image The cache matrix M ₁ of the initial value is zero, and its size is the same as the background image of the multispectral image, and the composite background image of the multispectral image is The gray value of each position of is directly assigned to the synthetic background image The corresponding position of cache matrix M ₁ ;

7.2) Statistically synthesized background images The maximum and minimum values of all values in the cache matrix M ₁ are denoted as and

7.3) According to the maximum and minimum Computationally synthesized background image Threshold T ₁ :

${\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; m</span>}}_{{\mathrm{<span\; class="notranslate">\; min</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 0.5</span>}<span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; m</span>}}_{{\mathrm{<span\; class="notranslate">\; max</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; m</span>}}_{{\mathrm{<span\; class="notranslate">\; min</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ;</span>$

7.4) Synthesize the background image of the multispectral image Compare the gray value of each pixel with its threshold T ₁ , classify the pixels whose gray value is greater than its threshold T ₁ into one category, and record the area composed of such pixels as S ₁₁ ; _The pixels of the threshold T1 are classified into another category, and the area composed of such pixels is denoted as _S12 ;

7.5) Calculate the mean value of all pixel gray values in regions S ₁₁ and S ₁₂ respectively, denoted as and Calculate the threshold R _{1 of the compressed value range of the cache matrix M 1 :}

${\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; A</span>}}_{{\mathrm{<span\; class="notranslate">\; S</span>}}_{\mathrm{<span\; class="notranslate">\; 11</span>}}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; A</span>}}_{{\mathrm{<span\; class="notranslate">\; S</span>}}_{\mathrm{<span\; class="notranslate">\; 12</span>}}}}{{\mathrm{<span\; class="notranslate">\; m</span>}}_{{\mathrm{<span\; class="notranslate">\; max</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; m</span>}}_{{\mathrm{<span\; class="notranslate">\; min</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; B</span>}<span\; class="notranslate">\; )</span>$

7.6) Normalize all values in the cache matrix M ₁ to obtain a normalized matrix The normalization matrix All values in are in the [0,1] interval;

7.7) For normalized matrix Each value in is multiplied by the threshold R ₁ to compress the range of values, so that the normalized matrix All the values in are compressed into the range [0, R ₁ ], and the compressed matrix is the weight matrix w ₁ .

Step 8: Establish a second weight matrix w ₂ according to the background image B ₂ of the panchromatic image.

8.1) Establish the buffer matrix M ₂ of the background image B ₂ of the full-color image, the initial value is zero, and its size is the same as the background image of the full-color image, and directly assign the gray value of each position of the background image B ₂ Give it the corresponding position of cache matrix M ₂ ;

8.2) The maximum and minimum values of all the values in the cache matrix M ₂ of the background image B ₂ of the panchromatic image are statistically denoted as and

8.3) According to the maximum and minimum Calculate the threshold T ₂ of the background image B ₂ of the panchromatic image,

${\mathrm{<span\; class="notranslate">\; T</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; m</span>}}_{{\mathrm{<span\; class="notranslate">\; min</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 0.5</span>}<span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; m</span>}}_{{\mathrm{<span\; class="notranslate">\; max</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; m</span>}}_{{\mathrm{<span\; class="notranslate">\; min</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ;</span>$

8.4) Compare the gray value of each pixel in the background image B ₂ of the full-color image with its threshold T ₂ , classify the pixels whose gray value is greater than the threshold T ₂ into one category, and record such pixels as regions Make S ₂₁ ; classify pixels whose gray value is smaller than the threshold T ₂ into another category, and record the area formed by these pixels as S ₂₂ ;

8.5) Calculate the mean value of the gray value of all pixels in the regions S ₂₁ and S ₂₂ respectively, denoted as and Calculate the threshold R _{2 of the compressed value range of the cache matrix M 2 :}

${\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; A</span>}}_{{\mathrm{<span\; class="notranslate">\; S</span>}}_{\mathrm{<span\; class="notranslate">\; twenty\; one</span>}}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; A</span>}}_{{\mathrm{<span\; class="notranslate">\; S</span>}}_{\mathrm{<span\; class="notranslate">\; twenty\; two</span>}}}}{{\mathrm{<span\; class="notranslate">\; m</span>}}_{{\mathrm{<span\; class="notranslate">\; max</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; m</span>}}_{{\mathrm{<span\; class="notranslate">\; min</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}}}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; B</span>}<span\; class="notranslate">\; )</span>$

8.6 ₎ Normalize all the values in the cache matrix _M2 of the background image B2 of the panchromatic image to obtain a normalized matrix The normalization matrix All values in are in the [0,1] interval;

8.7) For normalized matrix Each value in is multiplied by the threshold R ₂ to compress the value range, so that the normalized matrix All the values in are compressed to the range of [0, R ₂ ], and the compressed matrix is the weight matrix w ₂ .

Step 9: Multiply each direction subband coefficient z ₁ , z ₂ , ..., z _m of the panchromatic image I ₂ by the weight matrix w ₁ and w ₂ to obtain the direction of the fused first component image The sub-band coefficients l ₁ , l ₂ , ..., l _m directly assign the low-frequency coefficient x ₁ of the first principal component image to the low-frequency coefficient k of the fused first principal component image.

Step 10, performing an inverse Shearlet transform on the low-frequency coefficient k and multiple directional subband coefficients l ₁ , l ₂ , ..., l _m of the fused first principal component image to obtain the fused first principal component image

Step 11, the fused first principal component image and other component images obtained in step 4 except the first principal component Form a new data set, and perform inverse PCA transformation on the new data set to obtain the fused image I, as shown in Figure 3, where Figure 3(a) is the fused image enhanced by the weight matrix, and Figure 3 ( b) is a partially enlarged view of Fig. 3(a).

The above description is a specific example of the present invention. Obviously, for those skilled in the art, after understanding the content and principle of the present invention, it is possible to make changes in form and details without departing from the principle and structure of the present invention. Amendments and changes, but these amendments and changes based on the idea of the present invention are still within the protection scope of the claims of the present invention.

Claims (2)

1., with the multispectral of thin cloud and a panchromatic image fusion method, comprise the steps:

1) carry out up-sampling to the multispectral image with thin cloud, the size of the multispectral image after sampling is identical with the size of the full-colour image with thin cloud;

2) respectively down-sampling, Gassian low-pass filter and up-sampling are carried out successively to the multispectral and full-colour image after sampling, obtain the background image B of multispectral image ₁with the background image B of full-colour image ₂, wherein the wave band number of the background image of multispectral image is identical with the wave band number of multispectral image;

3) operation of thin cloud is removed respectively to each wave band of multispectral image and full-colour image:

3.1) calculate the gray average of multispectral image each band image and the gray average of full-colour image, namely the gray-scale value of pixels all in image is sued for peace, then divided by the sum of all pixels of image;

3.2) to each wave band of multispectral image, subtracting background image B ₁corresponding wave band, add the gray average of corresponding wave band, obtain the multispectral image I removing thin cloud ₁;

3.3) to full-colour image, its background image B is deducted ₂, add the average of its gray scale, obtain the full-colour image I removing thin cloud ₂;

4) to the multispectral image I removed after thin cloud ₁carry out PCA conversion, obtain each component image wherein represent the i-th principal component that multispectral image obtains after PCA conversion, i=1,2 ..., n, n are the sum of component;

5) to the first factor image obtained after PCA conversion with the full-colour image I removing thin cloud ₂, carry out Shearlet conversion respectively and decompose, by the first factor image be decomposed into a low frequency coefficient x ₁with multiple directions sub-band coefficients y ₁, y ₂..., y _m, will the full-colour image I of thin cloud be removed ₂be decomposed into a low frequency coefficient x ₂with multiple directions sub-band coefficients z ₁, z ₂..., z _m;

6) to step 2) background image of each wave band of multispectral image that obtains, calculate its average gray by following formula, obtain the synthesis background image B ' of multispectral image ₁:

Wherein represent background image B ₁the gray-scale value of kth wave band at coordinate (p, q) place, B ' ₁(p, q) represents the background image B ' of synthesis ₁at the gray-scale value at coordinate (p, q) place, N is the wave band number of multispectral image;

7) according to the synthesis background image B ' of multispectral image ₁, set up weight matrix w ₁:

7.1) synthesis background image B ' is set up ₁buffer memory matrix M ₁, initial value is zero, and its size is identical with the background image of multispectral image, by the synthesis background image B ' of multispectral image ₁the gray-scale value of each position, indirect assignment gives its buffer memory matrix M ₁relevant position;

7.2) statistics synthesis background image B ' ₁buffer memory matrix M ₁the maximal value of middle all values and minimum value, be denoted as respectively with

7.3) according to maximal value and minimum value calculate synthesis background image B ' ₁threshold value T ₁:

7.4) by the synthesis background image B ' of multispectral image ₁in the gray-scale value of each pixel and its threshold value T ₁compare, gray-scale value is greater than its threshold value T ₁pixel be divided into a class, and the region that this kind of pixel forms is denoted as S ₁₁; Gray-scale value is less than its threshold value T ₁pixel be divided into another kind of, and the region that this kind of pixel forms is denoted as S ₁₂;

7.5) difference zoning S ₁₁and S ₁₂in the average of all grey scale pixel values, be denoted as respectively with calculate synthesis background image B ' ₁buffer memory matrix M ₁the threshold value R of compression span ₁:

；

7.6) to described buffer memory matrix M ₁in all values be normalized, obtain normalization matrix M ' ₁, this normalization matrix M ' ₁in all value all in [0,1] interval;

7.7) to normalization matrix M ' ₁in each value be multiplied by its compression span threshold value R ₁, make normalization matrix M ' ₁in all values be compressed to [0, R ₁] in scope, the matrix after compression is exactly weight matrix w ₁;

8) according to the background image B of full-colour image ₂, set up weight matrix w ₂:

8.1) the background image B of full-colour image is set up ₂buffer memory matrix M ₂, initial value is zero, and its size is identical with the background image of full-colour image, by background image B ₂the gray-scale value of each position, indirect assignment gives its buffer memory matrix M ₂relevant position;

8.2) the background image B of full-colour image is added up ₂buffer memory matrix M ₂the maximal value of middle all values and minimum value, be denoted as respectively with

8.3) according to maximal value and minimum value calculate the background image B of full-colour image ₂threshold value T ₂,

8.4) by the background image B of full-colour image ₂in the gray-scale value of each pixel and its threshold value T ₂compare, gray-scale value is greater than its threshold value T ₂pixel be divided into a class, and this kind of pixel compositing area is denoted as S ₂₁; Gray-scale value is less than its threshold value T ₂pixel be divided into another kind of, and the region that these pixels form is denoted as S ₂₂;

8.5) difference zoning S ₂₁and S ₂₂in the average of all grey scale pixel values, be denoted as respectively with calculate the background image B of full-colour image ₂buffer memory matrix M ₂the threshold value R of compression span ₂:

；

8.6) to described buffer memory matrix M ₂in all values be normalized, obtain normalization matrix M ' ₂, this normalization matrix M ' ₂in all value all in [0,1] interval;

8.7) to normalization matrix M ' ₂in each value be multiplied by its compression span threshold value R ₂, make normalization matrix M ' ₂in all values be compressed to [0, R ₂] in scope, the matrix after compression is exactly weight matrix w ₂;

9) to full-colour image I ₂each directional subband coefficient z ₁, z ₂..., z _m, be multiplied by weight matrix w ₁and w ₂, obtain the directional subband coefficient l of the first component image merged ₁, l ₂..., l _m, by the low frequency coefficient x of the first factor image ₁as the low frequency coefficient k merging rear the first factor image;

10) to low frequency coefficient k and the multiple directions sub-band coefficients l of the first factor image after fusion ₁, l ₂..., l _m, carry out inverse Shearlet conversion, obtain the first factor image after merging

11) by the first factor image after fusion with step 4) other component images except the first factor of obtaining form new data set, and inverse PCA conversion is carried out to this new data set, obtain the image I after merging.

2. according to claim 1 with the multispectral of thin cloud and panchromatic image fusion method, wherein step 5) described in the first factor image with the full-colour image I removing thin cloud ₂carry out Shearlet decomposition respectively, carry out as follows:

5.1) carrying out Scale Decomposition by non-lower sampling Laplacian Pyramid Transform to image, is a low frequency coefficient and multiple high frequency coefficient by picture breakdown, and the number Sum decomposition number of plies of high frequency coefficient is identical, and determination Decomposition order of the present invention is 4 layers;

5.2) by Shear wave filter, each high frequency coefficient is decomposed into multiple directions, the present invention determines respectively from thick yardstick to thin yardstick, each high frequency coefficient is decomposed into 6,6,10 and 10 directions, obtains multiple directions sub-band coefficients.