An Improved IHS and Wavelet Transform Fusion Method for Nighttime Remote Sensing Images

Abstract

The present invention belongs to the technical field of remote sensing image fusion, and particularly relates to a method for fusing night-time remote sensing images by improving IHS and wavelet transform. The present invention fuses the "quasi-NPP-VIIRS" night-time light dataset and Landsat multispectral image data, performs multi-band synthesis on the night-time light remote sensing data, improves the image resolution (within 50 m), and realizes color night-time light remote sensing images. The present invention obtains the multi-band fused image after fusion, obtains the I component of the image brightness, fuses it with the corresponding panchromatic image to obtain a new I component of the image, and performs an IHS inverse transform on the H component and the S component to obtain the fused image. The I component contains details such as spatial structure and features. By using the present invention, detail information can be retained to a greater extent, and the accuracy of the fused image is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of remote sensing image fusion, and particularly relates to a method for fusing night light remote sensing images by improving IHS and wavelet transform. Background Technique

[0002] Traditional night light data is black and white images. The mainstream DMSP / OLS and NPP / VIIRS night light data have the deficiency of low resolution. By using multi-source image fusion algorithms, color night light remote sensing images can be obtained, which can further provide scientific references and services for fields such as urban light pollution and urbanization expansion.

[0003] At present, the multi-band fusion of night light remote sensing is still in the experimental stage, and multi-source remote sensing image fusion involves multiple fields. The current multi-source remote sensing image fusion is a fusion algorithm designed based on actual problem requirements and applications. There are situations such as spectral distortion, spatial information distortion, and high complexity of wavelet transform algorithms in traditional image fusion algorithms; with the rapid development of the application of remote sensing technology, existing research results still cannot meet actual needs, especially when effectively fusing multi-source remote sensing images between different sensors. There is still room for further expansion of multi-source remote sensing image fusion. In addition, each remote sensing image has certain limitations. Without changing the existing climate, make full use of the complementary information of different remote sensing images, and improve the spatio-temporal resolution, spectral resolution and other characteristics of the image according to the fusion algorithm. In recent years, with the improvement of computer performance, hardware technology and wavelet transform, wavelet image fusion algorithms have become a hot topic.

[0004] Although the traditional IHS transform algorithm for image fusion tries to retain as much information such as the hue and saturation of the multi-spectral image and the spatial detail features of the panchromatic image as possible, there are generally problems such as spectral distortion and spatial information distortion in the fusion process; the traditional wavelet transform decomposes the image into different frequencies for selective fusion, retaining the spectral and spatial structure features of the image to a great extent, but there are still problems with high algorithm complexity such as the number of wavelet decomposition levels and the selection of wavelet bases.

[0005] In the prior art, the publication number is CN114331936A, and the name is a method for fusing remote sensing images based on wavelet decomposition and improved IHS algorithm, which combines the improved HIS and wavelet transform. However, its I component undergoes weighted average and histogram matching, ignoring many detail components and having low accuracy. And its wavelet algorithm has high complexity. Summary of the Invention

[0006] In order to solve the above problems, the present invention provides a method for fusing night light remote sensing images by improving IHS and wavelet transform. The specific technical solutions are as follows:

[0007] An improved method for fusing night-time remote sensing images using IHS and wavelet transform, comprising the following steps:

[0008] Step S1, collect data, including NPP / VIIRS night-time light remote sensing images, Landsat multispectral images and Landsat panchromatic images. The Landsat multispectral image is the MS image, and the Landsat panchromatic image is the PAN image;

[0009] Step S2, data preprocessing, perform reprojection and resampling on the collected data to ensure that their image sizes and dimensions are consistent;

[0010] Step S3, multi-band fusion, fuse the preprocessed NPP / VIIRS night-time light remote sensing image with the preprocessed MS image to generate a color night-time light remote sensing image, i.e., a multi-band fusion image;

[0011] Step S4, IHS transformation of the image, perform IHS transformation on the multi-band fusion image processed in Step S3, transform the image color space to HIS, that is, decompose the original color image into three channels of R, G, and B, and extract the I component containing detail information and the H and S components containing spectral information of the multi-band fusion image in the IHS color space according to the conversion relationship between the IHS model and the RGB model;

[0012] Step S5, histogram matching, keep the H and S components of the multi-band fusion image unchanged, and perform histogram matching on the I component of the multi-band fusion image with spatial details and the preprocessed PAN image to obtain the PAN new image and I new component;

[0013] Step S6, two-dimensional wavelet decomposition and inverse transformation, decompose and perform inverse transformation on the I new component and the PAN new image separately by two-dimensional wavelet decomposition to obtain the I new ' component;

[0014] Step S7, image fusion, perform IHS inverse transformation on the obtained I new ' component, the H component, and the S component of the multi-band fusion image obtained in Step S4 to obtain a fused image.

[0015] Preferably, in Step S5, it specifically includes: calculating the gray value of the PAN image, performing histogram equalization on the PAN image and the I component of the multi-band fusion image, as shown in the following formula (1), and adjusting the gray level of the PAN image according to the corresponding relationship between Z k , P z to obtain an I new component with a higher matching degree to the original image; Inew The components are matched with the histogram of the PAN image to obtain the PAN new image;

[0016]

[0017] In the formula, Z k represents the gray value of the PAN image, and P(Z j ) represents the estimated value of the image gray probability when k = j. S r (Z k ) represents the sum of P(Z j ); L - 1 is the degree of freedom.

[0018] Preferably, the step S6 specifically includes the following steps:

[0019] Step S61: The I new component is decomposed by rows to obtain the low-frequency component L and the high-frequency component H of the I new component, and then decomposed by columns to obtain the low-frequency sub-component LL new of the I I , the detail feature LH I in the horizontal direction, the detail feature HL I in the vertical direction, and the diagonal feature HH I ;

[0020] Step S62: The PAN new image is decomposed by rows to obtain the low-frequency component L and the high-frequency component H of the I new component, and then decomposed by columns to obtain the low-frequency sub-component LL new of the PAN P image, the detail feature LH P in the horizontal direction, the detail feature HL P in the vertical direction, and the diagonal feature HH P ;

[0021] Step S63: Replace the low-frequency sub-component LL new of the PAN P image with the low-frequency sub-component LL new of the I I component to form a new PAN new image, and decompose the low-frequency sub-component LL new of the new PAN P image by rows to obtain the low-frequency component L and the high-frequency component H of the low-frequency sub-component LL new of the new PAN P image, and then decompose by columns to obtain the low-frequency sub-component LLLL new of the low-frequency sub-component LL P of the new PAN P, the detailed feature LLLH in the horizontal direction P , the detailed feature LLHL in the vertical direction P , the diagonal feature LLHH P ;

[0022] Step S64, perform two-dimensional inverse wavelet transform on the low-frequency sub-component LL of the new PAN new image, and finally obtain the I P ' component. new ' component.

[0023] Preferably, in the two-dimensional inverse wavelet transform of step S64, different fusion coefficients are assigned to the high-frequency component and the low-frequency component, where the high-frequency component adopts the absolute value strategy fusion coefficient and the low-frequency component adopts the average value strategy fusion coefficient.

[0024] Preferably, the absolute value strategy fusion coefficient for the high-frequency component is specifically shown in formula (2), and the average value strategy fusion coefficient for the low-frequency component is shown in formula (3):

[0025]

[0026]

[0027] where P = (m, n) represents the element value at the subscript (m, n) in the coefficient matrix;

[0028] C(F, p) represents the coefficient matrix of the image F in the P neighborhood space, and W max represents the maximum weight value, and W min represents the minimum weight value, and p = (m, n) represents the spatial position of the coefficient matrix.

[0029] Preferably, the wavelet basis function in step S6 is selected from the Haar wavelet function in the two-dimensional discrete wavelet transform.

[0030] Preferably, in step S2, the collected NPP / VIIRS nighttime light remote sensing image is resampled to the same spatial resolution as the Landsat multispectral image and / or the Landsat panchromatic image by using the cubic convolution interpolation method.

[0031] Preferably, step S3 further includes performing radiometric calibration and atmospheric correction processing on the multi-band fusion image.

[0032] The beneficial effects of the present invention are as follows: The present invention integrates the "similar NPP-VIIRS" nighttime light dataset and Landsat multispectral image data, performs multi-band synthesis on the nighttime light remote sensing data, improves the image resolution (within 50m), and realizes a color nighttime light remote sensing image. The present invention obtains the multi-band fused image after fusion, obtains the I component of the image brightness, and fuses it with the corresponding panchromatic image to obtain a new I component of the image, and performs an IHS inverse transformation on the H component and the S component to obtain the fused image. The I component contains detailed information such as spatial structure and features. By using the present invention, the detailed information can be retained to a greater extent, and the accuracy of the fused image is improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0033] In order to more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the following will briefly introduce the drawings required for the description of the specific embodiments or the prior art. In all the drawings, similar elements or parts are generally denoted by similar reference numerals. In the drawings, the elements or parts are not necessarily drawn to actual scale.

[0034] Figure 1 It is a schematic flowchart of the present invention;

[0035] Figure 2 It is a specific flowchart of the present invention;

[0036] Figure 3 It is a flowchart of the inverse wavelet transform of the present invention; DETAILED DESCRIPTION OF THE EMBODIMENTS

[0037] The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the drawings in the embodiments of the present invention. Obviously, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0038] It should be understood that when used in this specification and the appended claims, the terms "comprises" and "comprising" indicate the presence of the described features, wholes, steps, operations, elements, and / or components, but do not preclude the presence or addition of one or more other features, wholes, steps, operations, elements, components, and / or their combinations.

[0039] It should also be understood that the terms used in the specification of the present invention are only for the purpose of describing specific embodiments and are not intended to limit the present invention. As used in the specification of the present invention and the appended claims, unless the context clearly indicates otherwise, the singular forms "a", "an", and "the" are intended to include the plural forms.

[0040] It should also be further understood that the term "and / or" used in the specification and appended claims of the present invention refers to any combination and all possible combinations of one or more of the associated listed items, and includes these combinations.

[0041] As Figure 1-2 shown, the specific implementation manner of the present invention provides an improved method for fusing night-time remote sensing images of IHS and wavelet transform, including the following steps:

[0042] Step S1, collect data, including class NPP / VIIRS night-time light remote sensing images, Landsat multispectral images and Landsat panchromatic images. The Landsat multispectral image is the MS image, and the Landsat panchromatic image is the PAN image.

[0043] Step S2, data preprocessing, perform reprojection and resampling on the collected data to ensure that their image sizes and dimensions are consistent; among them, for the collected class NPP / VIIRS night-time light remote sensing images, use the cubic convolution interpolation method to resample the class NPP / VIIRS night-time light remote sensing images to the same spatial resolution as the Landsat multispectral images and Landsat panchromatic images.

[0044] Step S3, multi-band fusion, extract the near-infrared band and red band (Landsat7 ETM+ Band3-4, Landsat 8 OLI_TIRS Band4-5) of the preprocessed MS image and fuse them with the single-band class NPP / VIIRS night-time light remote sensing image to generate a color night-time light remote sensing image, that is, a multi-band fusion image; it also includes performing radiometric calibration and atmospheric correction processing on the multi-band fusion image.

[0045] Radiometric calibration refers to the process of correcting the random radiation distortion or aberration caused by external factors and correcting the image aberration. The radiometric calibration of Landsat data is to correct the DN value of the image and convert the original DN value into radiance value to eliminate errors such as sensor performance and solar altitude angle. Its radiance L λ The calculation method is expressed as:

[0046] L λ = M L × Q cal + A L ; (1)

[0047] In the formula: M L represents the image gain; A L represents the image offset.

[0048] Atmospheric correction is the process of eliminating the attenuation effect and backscattering caused by atmospheric scattering, resulting in differences in irradiance values and retrieving true physical model parameters such as surface emissivity. The FLAASH model is used for data correction, and the formula is shown as formula (2) below:

[0049]

[0050] In the formula: L represents the radiance of each pixel received, ρ represents the surface reflectivity of each pixel, and ρ ε represents the average reflectivity of each pixel and its neighborhood, S represents the spherical albedo of the atmosphere, and L a represents the path radiance of the atmosphere, and B is a coefficient determined by both the surface underlying surface and the atmosphere.

[0051] When dust or haze appears in the atmosphere, the "adjacent pixel effect" will occur in formula (2), masking the true surface conditions, and significant interpolation will occur for ρ and ρ ε , and at this time, ρ ε can be obtained through formula (3):

[0052]

[0053] In the formula: L ε is the average reflectivity of the specified pixel and its neighborhood.

[0054] Step S4, IHS transformation of the image. Perform IHS transformation on the multi-band fusion image processed in step S3, transform the image color space to IHS, that is, decompose the original color image into three channels of R, G, and B, and extract the I component containing detailed information and the H and S components containing spectral information of the multi-band fusion image in the IHS color space according to the conversion relationship between the IHS model and the RGB model.

[0055] Step S5, histogram matching. Keep the H and S components of the multi-band fusion image unchanged, and perform histogram matching on the I component of the multi-band fusion image with spatial details and the pre-processed PAN image to obtain the PAN new image and the I new component; specifically including: calculating the gray value of the PAN image, performing histogram equalization on the PAN image and the I component of the multi-band fusion image, as shown in formula (4) below, and adjusting the gray level of the PAN image according to the corresponding relationship between Z k , P z to obtain the I new component with a higher matching degree to the original image; the I new component is histogram-matched with the PAN image to obtain the PAN new image;

[0056]

[0057] In the formula, Z k represents the gray value of the PAN image, and P(Z j ) represents the estimated value of the image gray probability when k = j. S r (Z k ) represents the sum of P(Z j ); L - 1 is the degree of freedom.

[0058] Step S6: Two - dimensional wavelet decomposition and inverse transformation. Decompose and perform inverse transformation on the I new component obtained in step S4 and the PAN new image to obtain the I new ' component. As Figure 3 shown, it specifically includes the following steps:

[0059] Step S61: Perform row decomposition on the I new component to obtain the low - frequency component L and high - frequency component H of the I new component, and then perform column decomposition to obtain the low - frequency sub - component LL new of the I I , the detail feature LH I in the horizontal direction, the detail feature HL I in the vertical direction, and the diagonal feature HH I ;

[0060] Step S62: Perform row decomposition on the PAN new image to obtain the low - frequency component L and high - frequency component H of the I new component, and then perform column decomposition to obtain the low - frequency sub - component LL new of the PAN P image, the detail feature LH P in the horizontal direction, the detail feature HL P in the vertical direction, and the diagonal feature HH P ;

[0061] Step S63: Replace the low - frequency sub - component LL new of the PAN P image with the low - frequency sub - component LL new of the I I component to largely preserve the spectral information and spatial detail information of the image, form a new PAN new image, and perform row decomposition on the low - frequency sub - component LL new of the new PAN P image to obtain the low - frequency component L and high - frequency component H of the low - frequency sub - component LL new of the new PAN P image, and then perform column decomposition to obtain the low - frequency sub - component LL new of the new PAN PThe low-frequency subcomponent LLLL P , horizontal detail features LLLH P , vertical detail features LLHL P , diagonal feature LLHH P ;

[0062] Step S64, the new PAN new The low-frequency subcomponent LL of the image P Perform two-dimensional wavelet inverse transform and finally get I new ' component. The wavelet basis function is selected from the Haar wavelet function in the two-dimensional discrete wavelet transform.

[0063] In the two-dimensional wavelet inverse transform, different fusion coefficients are assigned to high-frequency components and low-frequency components, where the absolute value strategy is used for the high-frequency components and the average value strategy is used for the low-frequency components.

[0064] The absolute value strategy fusion coefficient of the high-frequency component is shown in formula (5), and the average value strategy fusion coefficient of the low-frequency component is shown in formula (6):

[0065]

[0066]

[0067] Where P = (m, n) represents the element value of the coefficient matrix with the table below (m, n);

[0068] C(F,p) represents the spatial coefficient matrix of image F in the P domain, W max represents the maximum weight, W min represents the minimum weight, and p=(m,n) represents the spatial position of the coefficient matrix.

[0069] Step S7, image fusion, the obtained I new The ' component is subjected to an IHS inverse transformation with the H component and the S component of the multi-band fused image obtained in step S4 to obtain a fused image.

[0070] Those of ordinary skill in the art will appreciate that the units of each example described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, computer software, or a combination of the two. In order to clearly illustrate the interchangeability of hardware and software, the composition of each example has been generally described in terms of function in the above description. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professional and technical personnel can use different methods to implement the described functions for each specific application, but such implementation should not be considered to be beyond the scope of the present invention.

[0071] In the embodiments provided in the present application, it should be understood that the division of units is only a logical function division. In actual implementation, there may be other division methods. For example, multiple units can be combined into one unit, one unit can be split into multiple units, or some features can be ignored, etc.

[0072] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements on some or all of the technical features. These modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention, and they should all be covered by the scope of the claims and the description of the present invention.

Claims

1. An improved IHS and wavelet transform-based fusion method for night light remote sensing images, characterized in that Including the following steps: Step S1, data collection, including NPP / VIIRS nighttime light remote sensing images, Landsat multispectral images, and Landsat panchromatic images. The Landsat multispectral image is the MS image, and the Landsat panchromatic image is the PAN image; Step S2, data preprocessing, reprojecting and resampling the collected data to ensure that their image sizes and dimensions are consistent; Step S3, multi-band fusion, fusing the preprocessed NPP / VIIRS nighttime light remote sensing image with the preprocessed MS image to generate a color nighttime light remote sensing image, i.e., a multi-band fusion image; Step S4, IHS transformation of the image, performing IHS transformation on the multi-band fusion image processed in Step S3, transforming the image color space to IHS, that is, decomposing the original color image into three channels of R, G, and B, and extracting the I component containing detail information and the H and S components containing spectral information of the multi-band fusion image in the IHS color space according to the conversion relationship between the IHS model and the RGB model; Step S5, histogram matching. Keep the H and S components of the multi-band fused image unchanged, and perform histogram matching between the I component of the multi-band fused image with spatial details and the preprocessed PAN image to obtain the PAN image and component; Step S6, two-dimensional wavelet decomposition and inverse transformation, perform two-dimensional wavelet decomposition and inverse transformation on the component and the PAN image respectively to obtain components; Step S7, image fusion, perform IHS inverse transformation on the obtained component and the component of the multi-band fusion image obtained in step S4, to obtain a fused image; The specific steps in Step S6 include the following: Step S61, decompose the component by rows to obtain the low-frequency component L and high-frequency component H of the component, and then decompose by columns to obtain the I low-frequency sub-component LL I of the I component, the detail feature LH in the horizontal direction I , the detail feature HL in the vertical direction, and the diagonal feature HH; Step S62, perform row decomposition on the PAN image to obtain the low-frequency component L and high-frequency component H of the component, and then perform column decomposition to obtain the low-frequency sub-component LL of the PAN P image, the detail feature LH in the horizontal direction P the detail feature HL in the vertical direction P and the diagonal feature HH P ; Step S63, replace the low-frequency sub-component LL of the PAN image with the low-frequency sub-component LL of the P component to form a new PAN image. Then, perform row decomposition on the low-frequency sub-component LL of the new PAN I image to obtain the low-frequency component L and high-frequency component H of the new PAN image. Next, perform column decomposition on the low-frequency sub-component LL of the new PAN image to obtain the low-frequency sub-component LLLL P of the new PAN image, the horizontal detail feature LLH P , the vertical detail feature LLHL P P , and the diagonal feature LLHH P P P ;​​​​ Step S64, perform an inverse two-dimensional wavelet transform on the low-frequency sub-component LL of the new PAN image obtained in step S63, and finally obtain P a component; ​ In the inverse two-dimensional wavelet transform of Step S64, different fusion coefficients are assigned to the high-frequency components and the low-frequency components, where the high-frequency components use the absolute value strategy fusion coefficient and the low-frequency components use the average value strategy fusion coefficient; The high-frequency components using the absolute value strategy fusion coefficient are specifically shown in formula (2), and the low-frequency components using the average value strategy fusion coefficient are shown in formula (3): ；（2） ；（3） Among them, represents the element value with the subscripts (m, n) in the coefficient matrix; Represent an image In The domain space coefficient matrix, Represent the maximum weight value, Represent the minimum weight value, Represent the spatial position of the coefficient matrix.

2. An improved method for fusing night light remote sensing images by IHS and wavelet transform according to claim 1, characterized in that, Specifically included in the step S5 are: calculating the gray value of the PAN image, performing histogram equalization on the PAN image and the I component of the multi-band fusion image, as shown in the following formula (1), and according to , the corresponding relationship between them, adjusting the gray level of the PAN image to obtain a component with a higher matching degree to the original image; The component is matched with the PAN image histogram to obtain the PAN image; ；（1） Wherein, represents the gray value of the PAN image, represents the estimated value of the image gray probability when k = j, represents the sum of; is the degree of freedom.

3. An improved method for fusing night light remote sensing images by IHS and wavelet transform according to claim 1, characterized in that, In Step S6, the wavelet basis function is selected from the Haar wavelet function in the two-dimensional discrete wavelet transform.

4. An improved IHS and wavelet transform-based fusion method for night light remote sensing images according to claim 1, characterized in that, In Step S2, the cubic convolution interpolation method is used to resample the collected NPP / VIIRS nighttime light remote sensing image to the same spatial resolution as the Landsat multispectral image and the Landsat panchromatic image.

5. An improved IHS and wavelet transform-based fusion method for night-time remote sensing images according to claim 1, characterized in that Step S3 also includes radiometric calibration and atmospheric correction processing of the multi-band fusion image.

Images