CN104574346A - Optical remote sensing image decomposition algorithm - Google Patents

Abstract

The invention belongs to the field of remote sensing image processing, in particular to an optical remote sensing image decomposition algorithm. It includes: step 1: data reading; read the stored original image data, step 2: logarithmic transformation; take the logarithm of F(x, y), that is, f(x, y)=lg(F(x, y)), step 3: wavelet decomposition; perform discrete wavelet transform, the wavelet base is Haar wavelet, step 4: inverse transform; replace {H _p,t } with zero, and perform wavelet inverse transform with the low-frequency subgraph L _t to obtain For the low-frequency image L _i in the space domain of image f, step 5: fractal calculation; step 6: loop judgment; if D _i ≤ D _d , execute step 7, otherwise add 1 to the value of i, and then execute step 3 wavelet decomposition, step 7: Optimal scale judgment; Step 8: Image generation. The effect of the invention is that it overcomes the problem that the optimal decomposition scale cannot be accurately judged in the traditional image decomposition algorithm.

Description

A Decomposition Algorithm for Optical Remote Sensing Image

technical field

The invention belongs to the field of remote sensing image processing, in particular to an optical remote sensing image decomposition algorithm.

Background technique

At present, in the field of remote sensing, in order to obtain the characteristics of the spatial distribution of surface material components more accurately, the method of terrain correction is often used to eliminate the influence of terrain. For this reason, it has become a recognized idea to use wavelet multi-resolution analysis method to study the characteristics of targets of different scales on remote sensing images. Object classification and recognition with features at different scales. However, how to determine the level of decomposition is a difficult problem. The current image decomposition algorithm usually uses different levels of decomposition results, and is determined through experiments compared with actual features. This method is mainly empirical, not only time-consuming, but also lacks theoretical basis, so it is difficult to promote.

Contents of the invention

The purpose of the present invention is to provide an optical remote sensing image decomposition algorithm aiming at the defects of the prior art.

The present invention is achieved in that a kind of optical remote sensing image decomposition algorithm comprises the following steps:

Step 1: Data reading

Read the stored original image data, and the read image is the gray value, that is, the data obtained in this step is a three-dimensional array, in which two dimensions are the horizontal and vertical coordinates of the image, and the third dimension is the gray value corresponding to the coordinates. This step The obtained result is represented by (x, y, F(x, y)), where x, y are the horizontal and vertical coordinates of the image, and F(x, y) is the gray value of the point (x, y),

Step 2: Logarithmic transformation

Take the logarithm of F(x,y), that is, f(x,y)=lg(F(x,y)),

The result obtained in this step is represented by (x, y, f(x, y)), and in the subsequent representation, f is also used to represent the set of images formed by f(x, y).

Step 3: Wavelet decomposition

The gray value obtained in step 2 is subjected to discrete wavelet transform using the Mallat algorithm, and the wavelet base is the Haar wavelet.

In this step, f is decomposed to obtain wavelet transform coefficients {L _t , H _p,t }, L _t represents the low-frequency sub-image of image f at scale t, and H _p,t represents the high-frequency sub-image of image f at scale t in direction p , where p=1, 2, 3, p=1 represents the horizontal direction, p=2 represents the vertical direction, p=3 represents the diagonal direction,

Step 4: Inverse Transform

Substitute {H _p,t } with zero, perform wavelet inverse transformation with the low-frequency sub-image Lt, and obtain the low-frequency image L _i in the spatial domain of the image f. In the first calculation, i=t, and L _i represents the spatial domain at the scale i In the subsequent iterative calculation, the value of i is determined according to the subsequent steps. In the subsequent steps, the L _i image is used as the image to be judged L ₀ ,

Step 5: Fractal calculation

Use the Fraclab module of Matlab to read in the low-frequency image L ₀ and DEM, and use the box dimension tool to calculate. Under the premise of ensuring that the cross-correlation coefficient is equal to 1, take five consecutive points so that the maximum fitting error and the fitting point span ratio are the smallest The box count dimension calculated by the regression equation is the calculated box count dimension value, and the box count dimensions D _i and D _d of the low-frequency image L ₀ and DEM are obtained respectively,

Step 6: Loop Judgment

If D _i ≤ D _d , go to step 7, otherwise add 1 to the value of i, and then go to step 3 for wavelet decomposition,

Step 7: Judgment of the best scale

Judging the absolute value of D _i-1 -D _t and the absolute value of D _i -D _d , if Abs(D _i -Dd)≤Abs(D _i-1 -D _d ), the best decomposition scale is i, Otherwise, the best scale is i-1, and the best scale is assigned a value to I _best ,

The Abs means the absolute value,

Step Eight: Image Generation

Replace {H _p,t } with zero, t=1,2,...,I _best , p=1,2,3, perform wavelet inverse transformation with the low-frequency subgraph under the scale I _best , and then perform inverse logarithmic transformation to get Low-frequency images in the spatial domain - terrain submaps,

Use the high-frequency subgraph set {H _p,t }, t=1,2,...,I _best , p=1,2,3 and zero to perform wavelet inverse transformation, and then perform inverse logarithmic change to obtain the high frequency in the space domain image - lithology submap,

Output a high-frequency submap representing lithology and a low-frequency submap representing topography.

An optical remote sensing image decomposition algorithm as described above, wherein the DEM is a surface elevation image defined in the same spatial domain as the image f, which is also a three-dimensional array, and can be used (x, y, D(x, y) ), where x and y are the horizontal and vertical coordinates of the image respectively, and D(x, y) is the elevation value of the point (x, y).

The effect of using the present invention is: firstly, the present invention converts the information contained in the remote sensing image from the result of multiplicative operation to the result of additive operation by taking logarithm; and then uses wavelet function to decompose the image. In order to determine the decomposition level, the optimal wavelet transform scale is judged by calculating whether the fractal dimensions of the DEM and the low-frequency submap are equal (within a certain error range), and finally the image is decomposed into the low-frequency topographic submap and the low-frequency topographic submap in the spatial domain. The high-frequency lithology submap in the space domain overcomes the problem of inaccurate determination of the optimal decomposition scale in the traditional image decomposition algorithm. Therefore, the method of combining fractal dimension calculation with wavelet decomposition overcomes the shortcomings of needing to rely on experience to judge the scale of wavelet decomposition and the effect of terrain correction is not obvious, and can greatly improve the effect of image decomposition. Fusion has important significance and practical value. The idea of the decomposition algorithm can also be extended to other optical image decomposition applications.

Description of drawings

Fig. 1 is the basic flowchart of the inventive method;

Fig. 2 is the wavelet decomposition result of TM image in DS mining area.

Detailed ways

An optical remote sensing image decomposition algorithm, comprising the following steps:

Step 1: Data reading

Read the stored original image data, and the read image is the gray value, that is, the data obtained in this step is a three-dimensional array, in which two dimensions are the horizontal and vertical coordinates of the image, and the third dimension is the gray value corresponding to the coordinates. The result obtained in this step is represented by (x, y, F(x, y)), where x, y are the horizontal and vertical coordinates of the image, and F(x, y) is the gray value of the point (x, y).

Step 2: Logarithmic transformation

Take the logarithm of F(x,y), that is, f(x,y)=lg(F(x,y)).

In this step, the two factors that affect the remote sensing image - surface objects and terrain, are converted from multiplicative operations to additive operations.

The result obtained in this step is represented by (x, y, f(x, y)), and in the subsequent representation, f is also used to represent the set of images formed by f(x, y).

Step 3: Wavelet decomposition

The gray value obtained in step 2 is subjected to discrete wavelet transformation using the Mallat algorithm, and the wavelet base is Haar wavelet.

In this step, f is decomposed to obtain wavelet transform coefficients {L _t , H _p,t }, L _t represents the low-frequency sub-image of image f at scale t, and H _p,t represents the high-frequency sub-image of image f at scale t in direction p . Here p=1, 2, 3, p=1 represents the horizontal direction, p=2 represents the vertical direction, and p=3 represents the diagonal direction. The calculation results are shown in Figure 2.

Step 4: Inverse Transform

Replace {H _p,t } with zero, and perform wavelet inverse transformation with the low-frequency sub-image Lt to obtain the low-frequency image L _i in the spatial domain of the image f. In the first calculation, i=t, L _i represents the low-frequency image in the spatial domain at scale i, and the value of i is determined according to subsequent steps in subsequent iterative calculations. In this application, the L _i image is taken as the image to be judged L ₀ .

Step 5: Fractal calculation

Use the Fraclab module of Matlab to read in the low-frequency image L ₀ and DEM, and use the box dimension tool to calculate. Under the premise of ensuring that the cross-correlation coefficient is equal to 1, take five consecutive points so that the maximum fitting error and the fitting point span ratio are the smallest The box count dimension calculated by the regression equation is the calculated box count dimension value. The box count dimensions D _i and D _d of the low frequency image L ₀ and DEM are obtained respectively.

The DEM is a surface elevation image defined in the same spatial domain as the image f, which is also a three-dimensional array and can be represented by (x, y, D(x, y)), where x and y are the horizontal and vertical coordinates of the image respectively , D(x, y) is the elevation value of the point (x, y).

Step 6: Loop Judgment

If D _i ≤ D _d , execute step seven, otherwise add 1 to the value of i, and then execute step three for wavelet decomposition.

In a certain calculation, when decomposing to the fifth level, the fractal dimension of L ₀ is 2.16599, which is smaller than the fractal dimension of DEM 2.27822, so it stops at the fifth level.

Step 7: Judgment of the best scale

Judging the absolute value of D _i-1 -D _t and the absolute value of D _i -D _d , if Abs(D _i -Dd)≤Abs(D _i-1 -D _d ), the best decomposition scale is i, Otherwise the optimal scale is i-1. The best scale attaches value to I _best .

Said Abs means absolute value.

Step Eight: Image Generation

Replace {H _p,t } with zero, t=1,2,...,I _best , p=1,2,3. Perform wavelet inverse transformation with the low-frequency submap under the scale I _best , and then perform inverse logarithmic transformation to obtain the low-frequency image in the spatial domain—topographic submap.

Use the high-frequency subgraph set {H _p,t }, t=1,2,...,I _best , p=1,2,3 and zero to perform wavelet inverse transformation, and then perform inverse logarithmic change to obtain the high frequency in the space domain Image - Lithology submap.

Output a high-frequency submap representing lithology and a low-frequency submap representing topography.

Claims (2)

1. a remote sensing image decomposition algorithm, is characterized in that, comprises the steps:

Step one: digital independent

Read the raw image data stored, the image read out is gray-scale value, and the data that namely this step obtains are three-dimensional array, wherein bidimensional is the transverse and longitudinal coordinate of image, and the third dimension is the gray-scale value that coordinate is corresponding, the result (x that this step obtains, y, F(x, y)) represent, wherein x, y is respectively the transverse and longitudinal coordinate of image, F(x, y) be point (x, y) gray-scale value

Step 2: log-transformation

To F(x, y) take the logarithm, i.e. f (x, y)=lg (F (x, y)),

The result that this step obtains represents with (x, y, f (x, y)), also represents the set of the image that f (x, y) is formed in follow-up expression with f,

Step 3: wavelet decomposition

Carry out wavelet transform to the gray-scale value Mallat algorithm that step 2 obtains, wavelet basis is Haar wavelet transform,

This step is decomposed f, obtains wavelet conversion coefficient { L _t, H _p,t, L _trepresent the low frequency subgraph picture of yardstick t hypograph f, H _p,trepresent the high frequency subimage in p direction under yardstick t, p=1 here, 2,3, p=1 represents horizontal direction, and p=2 represents vertical direction, and p=3 represents angular direction,

Step 4: inverse transformation

{ H is replaced with zero _p,t, carry out wavelet inverse transformation with low frequency subgraph Lt, obtain the low-frequency image L of image f spatial domain _i, i=t, L when first time calculates _irepresent the low-frequency image of the spatial domain under i yardstick, in successive iterations calculates, the value of i is determined according to subsequent step, by L in subsequent step _iimage judges image L as waiting ₀,

Step 5: fractal calculation

Low-frequency image L is read in by the Fraclab module of Matlab ₀and DEM, and calculate with box meter dimension instrument, under ensureing that cross-correlation coefficient equals the prerequisite of 1, getting continuous 5 box meter dimensions that matching maximum error and match point span are obtained than regression equation calculation time minimum is calculated box meter dimension value, obtains low-frequency image L respectively ₀with the box meter dimension D of DEM _iand D _d,

Step 6: cycle criterion

If D _i≤ D _d, perform step 7, otherwise make i value add 1, then perform step 3 wavelet decomposition,

Step 7: best scale judges

Judge D _i-1-D _tabsolute value and D _i-D _dthe size of absolute value, if Abs (D _i-Dd)≤Abs (D _i-1-D _d), best decomposition scale is i, otherwise best scale is i-1, and best scale assignments is to I _best,

Described Abs represents absolute value,

Step 8: Computer image genration

{ H is replaced with zero _p,t, t=1,2 ..., I _best, p=1,2,3, with yardstick I _bestunder low frequency subgraph carry out wavelet inverse transformation, then carry out antilogarithm change and obtain low-frequency image---the landform subgraph of spatial domain,

With high frequency subgraph set { H _p,t, t=1,2 ..., I _best, p=1,2,3 and zero carry out wavelet inverse transformation, then carry out high frequency imaging---the lithology subgraph that antilogarithm change obtains spatial domain,

Export and represent the high frequency subgraph of lithology and represent the low frequency subgraph of landform.

2. a kind of remote sensing image decomposition algorithm as claimed in claim 1, it is characterized in that: described DEM is the earth's surface elevation map picture being defined in the same space territory with image f, is three-dimensional array equally, can (x be used, y, D(x, y)) represent, wherein x, y is respectively the transverse and longitudinal coordinate of image, D(x, y) be the height value of point (x, y).