CN106651847A - Remote sensing image fog information detection method - Google Patents

Abstract

The invention discloses a method for detecting cloud and fog information of a remote sensing image, which comprises the following steps: inputting a remote sensing image and obtaining coarse cloud and fog estimation information. The rough transmittance map of the image scene and the atmospheric skylight value of the scene are obtained by using the dark channel prior to obtain the fine transmittance map. According to the atmospheric multiple scattering image model, the atmospheric component in the model is extracted, and the component is superimposed with the rough cloud and fog estimation information to obtain the cloud and fog spatial intensity distribution map. According to the observation data obtained from the ground observation points and the values of the corresponding points in the cloud and fog spatial intensity distribution map, the logarithmic function is used for data fitting, so as to obtain the mapping relationship between the detected physical quantities and the cloud and fog spatial intensity distribution information required by the actual scene, and realize Quantitative measurement of cloud and fog atmospheric information. Based on the image processing, the present invention provides a convenient and effective detection method for urban air environment pollution monitoring by obtaining atmospheric cloud intensity and spatial distribution information of the scene from the image.

Description

A cloud and fog information detection method for remote sensing images

technical field

The invention relates to remote sensing detection technology and computer image processing calculation, in particular to a detection method for obtaining remote sensing image cloud information by using image processing.

Background technique

In today's society, smog and PM2.5 pollution are the key governance issues in most urban areas in China. However, various cities have adopted various emission reduction measures, including restricting cars with odd and even numbers, reducing carbon emissions, etc., but the overall air quality has not changed accordingly. The reason is that the monitoring and law enforcement in the process of air pollution control are not enough.

Usually, for the detection of air pollution such as smog, continuous observation of ground monitoring points is often used for monitoring. However, there are rich dynamic characteristics in the characteristics of atmospheric composition, and it is difficult to grasp the detailed situation and spatial distribution characteristics of air pollution in a large range only by relying on the data of a few monitoring points. At the same time, the atmospheric flow will also lead to the possibility that the continuous observation data may be polluted by the surrounding environment, making it difficult to determine the exact location of the air pollution source. Therefore, how to effectively grasp the location and intensity of air pollution sources and analyze the causes of smog in different cities is an important demand for current social environmental monitoring.

Contents of the invention

The technical problem to be solved by the present invention is to provide a technical means to obtain the corresponding high-resolution atmospheric cloud intensity and spatial distribution information for high-resolution satellite images of certain weather conditions.

In order to solve the above-mentioned technical problems, the technical solution adopted in the present invention is: a cloud and fog information detection technology for remote sensing images, and its implementation steps are as follows:

(1) Input a remote sensing image, and use low-pass filter and guided filter to obtain rough cloud and fog estimation map.

(2) Use the dark channel prior to obtain the rough transmittance map of the image scene and the atmospheric skylight value of the scene, based on the rough transmittance combined with the rough cloud estimation map estimated in step (1), use the guided filter to obtain the fine transmittance Overrate graph.

(3) Based on the filtered input image in step (1) and the fine transmittance and skylight values obtained in step (2), according to the atmospheric multiple scattering image model It is possible to calculate the restored image while also extracting the atmospheric component in the model This component is superimposed with the rough cloud and fog estimation information obtained in step (1), which is the cloud and fog spatial intensity distribution information extracted from the remote sensing image.

(4) According to the observation data obtained from the ground observation points in the remote sensing image scene and the value of the corresponding point in the cloud and fog spatial intensity distribution map obtained in step (3), use the logarithmic function to perform data fitting, so as to obtain the detection required for the actual scene The mapping relationship between physical quantities and the spatial intensity distribution information of clouds and fog realizes the quantitative measurement of cloud and fog atmospheric information.

Further, the detailed steps in the step (1) include:

a) Input a remote sensing image, perform Fourier transform on it, obtain its spectrogram, and construct the following low-pass filter to filter the image at the same time:

Where H(u, v) is the filter function, u, v are the frequency coordinates, and σ ₀ is the cut-off frequency. The non-uniform cloud background of the filtered output is as follows:

where B _cloud is the low-frequency information of the input image, I is the input image, with These are the Fourier transform and inverse Fourier transform operators.

b) Contrast expansion is performed on the preliminary estimated cloud and fog distribution information image, and the corresponding adjustment formula is as follows:

in is the threshold value of the adjustment formula, max(B _cloud ) and min(B _cloud ) are the maximum and minimum values of the initially estimated non-uniform cloud background image respectively, and k ₁ , k ₂ and λ are the parameters of the adjustment formula.

And carry out zero-frequency compensation on it to obtain the final and preliminary estimated atmospheric cloud distribution information. The compensation formula is as follows:

B' _cloud = B _cloud -offset (4)

Where offset is the compensation constant of zero-frequency compensation, and then use the input image as a guide map, and use the guide filter to process the obtained atmospheric cloud and fog distribution map to obtain a preliminary cloud and fog spatial distribution estimation map with edge preservation and smoothness.

c) Remove the preliminary estimated atmospheric cloud distribution information from the input image:

I'(x,y)=I(x,y)-B' _cloud (x,y) (5)

The obtained images are used to further optimize the extraction results of atmospheric cloud and fog distribution information in the next processing.

Further, the detailed steps in the step (2) include:

a) Extract the dark channel of the image finally obtained in step 1 according to the prior knowledge of the dark channel:

Among them, I _dark is the dark channel of the acquired image, I ^c represents each color channel component of the image, and Ω(x) represents the local neighborhood centered on the x pixel position.

b) Extract part of the brightest area pixels in the top 0.1% from the dark channel of the above image, and calculate the pixel mean value of each channel in this area as the skylight value A.

c) According to the nature of the dark channel, its value is close to 0 for the fog-free area. Therefore, the initial transmittance of the image is calculated as follows:

Where t' is the calculated initial transmittance, ω is a constant coefficient, and ^Ac is the channel component of the skylight value.

d) Use the preliminary cloud and fog spatial distribution estimation map calculated in step (1) as a guide map, and use the guide filter to fine-tune the obtained initial transmittance to obtain a fine and high-resolution scene transmittance map.

Further, according to the atmospheric multiple scattering model in the step (3), its form is as follows:

Among them, I is the image obtained by subtracting the preliminary cloud distribution estimation from the input image, and A and t are the skylight value and transmittance image obtained in step 2. APSF _o and APSF _a are the atmospheric point spread function and the skylight point spread function respectively, which are solved by using the generalized Gaussian distribution, and its form is as follows:

Where x, y is the image coordinate position, Γ(.) is the gamma function, p and σ are the atmospheric parameters and are calculated as follows:

p=kT(10)

In the formula, T is the atmospheric optical thickness, which is calculated according to the relationship between transmittance and atmospheric optical thickness t=e- ^T , k is a parameter constant, and q is the forward scattering factor, which is a constant related to weather conditions. The optical thickness value is -log(t) for the atmospheric point spread function, and -log(1-t) for the skylight point spread function. So we can calculate the atmospheric component in the model:

Finally, combined with the preliminary atmospheric distribution information obtained in step 1, we can obtain the final atmospheric cloud distribution information map:

B"' _cloud = B' _cloud + B" _cloud (13)

Further, in the step (4), according to the input remote sensing image and its scene imaging range, the actual atmospheric measurement data (such as PM2.5 concentration, air quality index, etc.) of the corresponding ground observation point during the imaging of the remote sensing image is obtained as sample data Y, then take the corresponding point in the cloud and mist spatial intensity distribution map obtained in step (3) as the sample data X, and use the exponential function of the following form to fit the above-mentioned obtained data:

y=-Alog(B(1-x))+C (14)

Among them, y is the atmospheric measurement data value, x is the obtained cloud and fog information detection map data value, and A, B, and C are the parameters that need to be fitted. After obtaining the fitting parameters, the image obtained in step (3) is mapped point-by-point, and finally the quantitative atmospheric measurement data distribution map is obtained.

The beneficial effects of the present invention are: the present invention proposes a method for detecting cloud and fog information of remote sensing images, which realizes obtaining cloud and fog distribution information of corresponding scenes by inputting remote sensing images under certain atmospheric conditions. In the present invention, only simple input of remote sensing images is required without introducing other measurement parameters, and the cloud and fog distribution information of the atmosphere can be obtained while making the remote sensing images clear through the algorithm of digital image processing.

Description of drawings

Fig. 1 is a schematic flow chart of an embodiment of the present invention.

Fig. 2 is the original image of the embodiment of the present invention.

Fig. 3 is a remote sensing image obtained by the embodiment of the present invention after removing the preliminary cloud and fog distribution.

FIG. 4 is the optimized image transmittance extracted and obtained by the embodiment of the present invention.

Fig. 5 is the result of spatial intensity distribution of atmospheric cloud and mist obtained by the embodiment of the present invention.

Fig. 6 is a diagram of assumed positions of ground observation stations according to an embodiment of the present invention.

Fig. 7 is a fitting curve with ground data points in an embodiment of the present invention.

detailed description

The present invention will be described in further detail below in conjunction with the accompanying drawings and specific embodiments.

As shown in Figure 1, the implementation steps of the remote sensing image defogging method under the non-uniform cloud and fog conditions in this embodiment are as follows:

(1) Input a pair of remote sensing images (see Figure 2), select the cut-off frequency parameter as σ ₀ =1.7, and construct a low-pass filter Perform a low-pass filter operation on the input image. Set the parameters k ₁ =5, k ₂ =5 (for 8-bit images) and λ=2 of the adjustment formula to adjust the output image of the low-pass filter:

Perform zero-frequency compensation on the image output by the adjustment function to obtain the final and preliminary estimated atmospheric cloud distribution information. The compensation formula is as follows:

B' _cloud = B _cloud -offset

Among them, offset is the compensation constant of zero-frequency compensation, and its value is the average value of each channel of the input image. Then, the input image is used as the guide map, and the obtained atmospheric cloud and fog distribution map is processed by the guide filter to obtain a preliminary cloud and fog space with edge preservation and smoothness. Distribution Estimation Plot. At the same time, the remote sensing image I'(x, y) that has been removed from the coarse-step cloud and fog distribution is obtained, and the result is shown in Figure 3.

(2) For the remote sensing image I'(x,y) obtained in step 1, the dark channel of the image is extracted according to the prior knowledge of the dark channel:

Among them, I _dark is the dark channel of the acquired image, I ^c represents the red, green and blue channel components of the image, and Ω(x) represents the local neighborhood centered on the x pixel position. Extract part of the brightest area pixels in the top 0.1% from the dark channel of the above image, and calculate the pixel average value of each channel in this area as the skylight value A, which is A ^R =255, A ^G =254, A ^B =250 in this embodiment.

According to the nature of the dark channel, its value is close to 0 for no fog area. Therefore, the initial transmittance of the image is calculated as follows:

Where t' is the calculated initial transmittance, ω is a constant coefficient, and ^Ac is the channel component of the skylight value.

Use the preliminary cloud and fog spatial distribution estimation map calculated in step (1) as a guide map, and use the guide filter to fine-tune the obtained initial transmittance to obtain a fine and high-resolution scene transmittance map, as shown in Figure 4.

(3) Then use the atmospheric multiple scattering model, its form is as follows:

Among them, I' is the image obtained by subtracting the preliminary cloud and fog distribution estimation from the input image, and A and t are the skylight value and transmittance image obtained in step 2. APSF _o and APSF _a are the atmospheric point spread function and the skylight point spread function respectively, which are solved by using the generalized Gaussian distribution, and its form is as follows:

Where x, y is the image coordinate position, Γ(.) is the gamma function, p and σ are the atmospheric parameters and are calculated as follows:

p=kT (17)

In the formula, T is the atmospheric optical thickness, which is calculated according to the relationship between transmittance and atmospheric optical thickness t=e- ^T , k is 0.5, and q is the forward scattering factor, and its value is shown in the following table. In this embodiment, is 0.7.

0.0-0.2 0.2-0.7 0.7-0.8 0.8-0.85 0.85-0.9 0.9-1.0 Air Aerosols Haze Misty Fog Rain

For the skylight point spread function, its optical thickness value is -log(1-t). So we can calculate the atmospheric component in the model:

Finally, combined with the preliminary atmospheric distribution information obtained in step 1, we can obtain the final atmospheric cloud distribution information map, as shown in Figure 5:

B"' _cloud = B' _cloud + B" _cloud (19)

(4) As shown in Figure 6, it is assumed that there are three ground atmospheric parameter observation points A, B, and C in the input remote sensing image area, and the air quality index (AQI) values detected during satellite imaging are 170, 142, and 155, respectively. The data of the corresponding point positions in the atmospheric cloud distribution information map obtained in step 3 are 0.4826, 0.3518, and 0.4132 respectively. Fit the above data using a logarithmic function of the form:

y=-Alog(B(1-x))+C (20)

The fitting parameters were A=124.0137, B=15.1922, C=425.8834. The fitting results are shown in Figure 7. Based on the fitted logarithmic function, we can calculate the air quality index at any position within the scene range of the remote sensing image. For example, the air quality index at D is 142 (the calculation result is rounded), and the result at E is 152.

Compared with the traditional method, this method does not require other data from satellite detection, and only realizes the cloud and fog layer extraction of the image through digital image processing technology, so as to achieve high-resolution cloud and fog spatial distribution and intensity map while making the remote sensing image clearer . At the same time, if equipped with ground observation point data, it can be further transformed from qualitative observation to quantitative observation.

Claims (5)

1. a cloud and fog information detection method of remote sensing image, it is characterized in that, the method comprises the following steps: (1) Input a remote sensing image, and use low-pass filter and guided filter to obtain rough cloud and fog estimation map. (2) Use the dark channel prior to obtain the rough transmittance map of the image scene and the atmospheric skylight value of the scene, based on the rough transmittance map combined with the rough cloud estimation map obtained in step 1, use the guided filter to obtain the fine transmittance picture. (3) Based on the input map processed in step 1 and the fine transmittance map and atmospheric skylight value obtained in step 2, according to the atmospheric multiple scattering image model While calculating the restored image, the atmospheric component in the model is extracted Superimpose this component with the rough cloud and fog estimation map obtained in step 1, which is the cloud and fog spatial intensity distribution map extracted from the remote sensing image. (4) According to the observation data obtained from the ground observation points in the remote sensing image scene and the value of the corresponding point in the cloud and fog spatial intensity distribution map obtained in step 3, use the logarithmic function to perform data fitting, so as to obtain the detection physical quantity and The mapping relationship between cloud and fog spatial intensity distribution information realizes the quantitative measurement of cloud and fog atmospheric information. 2. the cloud and fog information detection method of remote sensing image as claimed in claim 1, is characterized in that, described step 1 is specifically: Input a foggy remote sensing image, construct the following low-pass filter and adjustment function to preprocess the image: $\mathrm{<span\; class="notranslate">\; h</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; u</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; exp</span>}<span\; class="notranslate">\; \&lsqb;</span><span\; class="notranslate">\; -</span>\frac{{\mathrm{<span\; class="notranslate">\; u</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; v</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; \&sigma;</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; \&rsqb;</span>$ ${\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; u</span>}\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\left\{\begin{array}{c}{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; u</span>}\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{<span\; class="notranslate">\; \&lsqb;</span>\frac{{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; u</span>}\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}{\mathrm{<span\; class="notranslate">\; max</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; u</span>}\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}<span\; class="notranslate">\; \&rsqb;</span>}^{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}<span\; class="notranslate">\; \&times;</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; u</span>}\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ></span>{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}\\ {\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; u</span>}\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>{<span\; class="notranslate">\; \&lsqb;</span>\frac{{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; u</span>}\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span>}{{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; min</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; u</span>}\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; \&rsqb;</span>}^{\mathrm{<span\; class="notranslate">\; \&lambda;</span>}}<span\; class="notranslate">\; \&times;</span>{\mathrm{<span\; class="notranslate">\; k</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; c</span>}\mathrm{<span\; class="notranslate">\; l</span>}\mathrm{<span\; class="notranslate">\; o</span>}\mathrm{<span\; class="notranslate">\; u</span>}\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&le;</span>{\mathrm{<span\; class="notranslate">\; B</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}\end{array}\right.$ In order to obtain the rough fog distribution map of the image, and then use the input image combined with the guided filtering operation to perform edge-preserving and fine-grained operations on the fog rough distribution map to optimize the edge information of the cloud layer. 3. the cloud and fog information detection method of remote sensing image as claimed in claim 1, is characterized in that, described step 2 is specifically: Extract the dark channel of the image based on the dark channel prior knowledge: Use the dark channel to obtain the atmospheric skylight value A and image transmittance t of the image: Then use the rough cloud and fog estimation map obtained in step 1 as a guide map, and use the guided filtering operation to obtain the fine transmittance map of the scene. 4. the cloud and fog information detection method of remote sensing image as claimed in claim 1, is characterized in that, described step 3 is specifically: Based on atmospheric multiple scattering model Dehaze and restore the remote sensing image that has been filtered in step 1. Atmospheric components in the model are extracted during restoration And this component is superimposed with the rough cloud and fog estimation map extracted in step 1 to obtain the final cloud and fog spatial intensity distribution map of the remote sensing image. 5. the cloud and fog information detection method of remote sensing image as claimed in claim 1, is characterized in that, described step 4 is specifically: According to the actual atmospheric measurement data detected at the corresponding time at the corresponding ground observation point within the imaging range of the remote sensing image, combined with the corresponding points in the cloud and fog spatial intensity distribution map obtained in step 3, the exponential function of the following form is used for fitting: y=-Alog(B(1-x))+C Among them, y is the atmospheric measurement data value, x is the obtained cloud and fog spatial intensity distribution map data value, and A, B, and C are the parameters that need to be fitted. After obtaining the fitting parameters, the image obtained in step 3 is mapped point-by-point, and finally the quantitative atmospheric measurement data distribution map is obtained.