CN109300115B

CN109300115B - Object-oriented multispectral high-resolution remote sensing image change detection method

Info

Publication number: CN109300115B
Application number: CN201811018172.6A
Authority: CN
Inventors: 石爱业; 王维; 王鑫; 马贞立
Original assignee: Hohai University HHU
Current assignee: Hohai University HHU
Priority date: 2018-09-03
Filing date: 2018-09-03
Publication date: 2021-11-30
Anticipated expiration: 2038-09-03
Also published as: CN109300115A

Abstract

The invention discloses an object-oriented multispectral high-resolution remote sensing image change detection method, which comprises the following steps of: inputting two high-resolution optical remote sensing images in the same area and at different time phases; using remote sensing software ENVI to X₁And X₂Carrying out image registration, namely carrying out coarse correction and fine correction; method for detecting MAD by utilizing multivariate change to X₁And X₂Carrying out radiation normalization correction; segmenting based on SLIC algorithm; calculating a difference image of the multi-temporal image, clustering by using a K-Means algorithm, and determining a CRF unitary energy item according to a clustering result and in combination with a segmentation object of the SLIC; constructing a neighborhood system of the CRF according to the segmentation image of the SLIC, and constructing binary item energy of the CRF according to the difference image and the space coordinate of the object; and optimizing the CRF energy item U to obtain a final change detection result. The invention not only considers the spectrum difference of the adjacent objects, but also considers the space position relation of the adjacent objects, thereby improving the precision of change detection.

Description

Object-oriented multispectral high-resolution remote sensing image change detection method

Technical Field

The invention relates to an object-oriented multi-temporal multispectral high-resolution remote sensing image unsupervised change detection method based on a Conditional Random Field (CRF), and belongs to the technical field of optical remote sensing image change detection.

Background

With the continuous accumulation of multi-temporal remote sensing data and the successive establishment of spatial databases, how to extract and detect change information from the remote sensing data has become an important research topic of remote sensing science and geographic information science. According to the remote sensing images of the same region in different time phases, the information of dynamic changes of cities, environments and the like can be extracted, and scientific decision-making basis is provided for departments of resource management and planning, environmental protection and the like.

The change detection of the remote sensing image is to quantitatively analyze and determine the characteristics and the process of the surface change from the remote sensing data in different periods. Various scholars put forward a plurality of effective detection algorithms from different angles and application researches, and generally speaking, according to whether a training sample is needed in the detection process, the change detection can be divided into three main categories: unsupervised change detection algorithms, semi-supervised change detection algorithms and supervised change detection algorithms.

The pixel-based unsupervised change detection method is suitable for remote sensing images with medium-low resolution, and the methods imply that pixels have independence in space. For a high-resolution image, because the information such as the structure and texture is more prominent, the efficiency and the obtained result information are very limited by using the traditional processing mode based on the pixel spectral statistics, and the processing result often has many small plaques, so that the analysis and description of the small plaques are difficult. The object-oriented change detection is more beneficial to the combined utilization of knowledge, and can also more effectively utilize the multi-feature advantages of high-resolution images.

In the change detection modeling process, the pixel and the neighborhood thereof have high correlation, and the neighborhood relationship can be described by two probability map models, namely MRF and CRF. The CRF can simultaneously consider the spatial context information of the observation field and the marker field, is more flexible than MRF modeling, and has advantages in the field of multi-temporal phase change detection. For example, the SAR image change detection method based on non-stationary analysis and conditional random field (2015, chinese patent of invention, west ampere electronic technology university) proposed by lister et al, but the inventive technology is supervised type change detection, and a large amount of labor cost and the like are consumed in practical application to construct a training sample. The unsupervised transformation detection method based on CRF provided by Cao and the like is applied to the detection process of the multispectral multi-temporal remote sensing image, and the detection precision is improved. However, the method is based on pixels, and the efficiency and the accuracy of detection have certain limitations. Object-oriented chi-squared transform based detection methods (OBCT) proposed by descl é et al can overcome the limitations of pixel-based change detection, but do not take into account the correlation between objects.

Disclosure of Invention

Aiming at the defects in the prior art, the invention aims to provide an object-oriented multispectral high-resolution remote sensing image change detection method, which combines the pixel-based and object-oriented technologies to ensure that the change detection result is more stable, and not only the spectral difference of adjacent objects but also the spatial position relation of the adjacent objects are considered when CRF binary energy items are constructed, so that the change detection precision is improved.

In order to achieve the purpose, the invention is realized by the following technical scheme:

the invention discloses an object-oriented multispectral high-resolution remote sensing image change detection method, which comprises the following steps of:

step 1: inputting two high-resolution optical remote sensing images in the same area and different time phases, and respectively recording the images as: x₁And X₂；

Step 2: using remote sensing software ENVI to X₁And X₂Carrying out image registration, namely carrying out coarse correction and fine correction;

and step 3: method for detecting MAD by utilizing multivariate change to X₁And X₂Carrying out radiation normalization correction;

and 4, step 4: superposing the preprocessed multi-temporal high-resolution images, and then segmenting based on an SLIC algorithm;

and 5: calculating the difference image of the multi-temporal image, and recording the image composed of the difference images of each wave band as X_MThen, clustering by using a K-Means algorithm, and determining a CRF unitary energy item according to a clustering result;

step 6: constructing a multi-neighborhood system of the CRF according to the segmentation image of the SLIC, and constructing binary item energy of the CRF according to the difference image in the step 5;

and 7: and optimizing the CRF energy item U according to an optimization algorithm and a cyclic belief propagation (LBP) algorithm to obtain a final change detection result.

In step 2, for the coarse correction, the specific steps are as follows:

2.1) displaying the reference image and the image to be corrected;

2.2) collecting ground control points GCPs, wherein the GCPs are uniformly distributed in the whole image, and the number of the GCPs is at least more than or equal to 25;

2.3) calculating the root mean square error;

2.4) selecting a polynomial model;

2.5) resampling and outputting by using bilinear interpolation.

Obtaining the value of the unknown function f at the point P as (x, y), and obtaining the value of the unknown function f at Q based on the preset known function f₁₁＝(x₁,y₁),Q₁₂＝(x₁,y₂),Q₂₁＝(x₂,y₁) And Q₂₂＝(x₂,y₂) Values of four points, if a coordinate system is selected such that the coordinates of the four points are (0,0), (0,1), (1,0), and (1,1), respectively, then the bilinear interpolation formula is expressed as:

f(x,y)≈f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy(1)。

in step 2, the multispectral remote sensing image data which is to undergo geometric rough correction is subjected to geometric fine correction by utilizing an automatic matching and triangulation method.

The triangulation method is to construct a Delaunay triangulation network by adopting a point-by-point insertion method, determine affine transformation model parameters inside each triangle by utilizing the row and column numbers of three vertexes of each triangle and the geographic coordinates of the same-name points of the corresponding reference image, correct the image to be corrected and obtain the corrected remote sensing image.

In step 3, the specific method for the normalized correction of the MAD relative radiation is as follows:

firstly, finding out a linear combination of brightness values of each wave band of the two-stage image through typical correlation analysis to obtain a difference image with enhanced change information, and determining changed and unchanged areas by taking a threshold value for a normalized MAD variable square sum;

and then, completing relative radiation correction through a mapping equation of the two-time phase pixel corresponding to the unchanged area.

The step 4 is realized by the following steps:

firstly, overlapping the preprocessed two-time-phase multispectral high-resolution remote sensing images. Then the following steps are carried out:

4.1) initializing seed points, namely cluster centers: uniformly distributing seed points in the image according to the number of the preset super pixels; based on the total N pixel points of the preset picture, the picture is pre-divided into R super-pixels with the same size, the size of each super-pixel is N/R, and the distance between adjacent seed points, namely the step length, is

4.2) reselecting the seed point in the n multiplied by n neighborhood of the seed point;

4.3) distributing a class label to each pixel point in the neighborhood around each seed point, wherein the search range is limited to 2S multiplied by 2S;

4.4) for each searched pixel point, respectively calculating the color distance and the space distance between the pixel point and the seed point; taking the seed point corresponding to the minimum value as a clustering center of the pixel points;

4.5) iterating step 4.1) -step 4.4) until the error converges.

In step 4.2), the specific selection method of the seed points comprises the following steps:

and calculating gradient values of all pixel points in the neighborhood, and moving the seed point to the position with the minimum gradient in the neighborhood.

In step 5, based on the preset cluster number K being 2, the specific determination method for the CRF unitary energy term is as follows:

5.1) randomly selecting K clustering centers;

5.2) classifying each data point into the category represented by the central point closest to the data point;

5.3) recalculating centers of the obtained classes;

5.4) iterating the step 5.2) to the step 5.3) until the new center is equal to the original center or smaller than a set threshold value;

5.5) after the K-Means algorithm is converged, counting the number of objects which belong to different classes in each object, and recording the number of objects j which belong to non-variable classes as N_j1The number of the variants is N_j2Then the unary energy term of object jComprises the following steps:

in step 6, the binary term energy for constructing CRF is as follows:

in the formula, i, j represents two adjacent objects or nodes of the image (in the invention, each object is regarded as one node in the graph model); c. C_i,c_jRepresenting class labels at two nodes i, j; two categories are changed and unchanged; l (i, j) represents the length of the common boundary of the two objects; [ c ] is_i≠c_j]Represents a 0-1 indicator function; f_M(i)、F_M(j) Respectively representing vectors consisting of the mean of the difference image values and the mean of the coordinates inside the objects i and j, wherein F_M(i) The construction of (a) is as follows:

F_M(i)＝[Mean(Dif_i)Mean(A_ix)Mean(A_iy)] (4)

wherein Mean (Dif)_i) Representing the Mean of the magnitude of the difference image within the object i, Dif representing the Mean of the values of the difference image within the object i, Mean (A)_ix) Mean, Mean (A) representing the abscissa of all pixels within the object i_iy) Represents the mean of the ordinate of all pixels within the object i;

the final CRF energy term constructed was as follows:

wherein U represents the total energy term for CRF, U_j、U_i,jA unitary energy term representing an object j or a node j, respectively, a binary energy term between adjacent objects j and i, S represents a set of nodes, E represents a set of adjacent nodes, and parameter β is a regularization parameter for controlling the ratio between the unitary and binary energy termsAnd (4) heavy.

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

(1) when the CRF unitary energy item is constructed, firstly, the category of each pixel is obtained by utilizing the K-means clustering based on the pixel. Then, based on the object, the occurrence frequency of each pixel in the object is counted, a robust unary energy item is further constructed, and the pixel-based technology and the object-facing technology are combined, so that the change detection result is more robust.

(2) When the CRF binary energy item is constructed, not only the spectrum difference of adjacent objects but also the spatial position relation of the adjacent objects are considered, so that the accuracy of change detection is improved.

Drawings

FIG. 1 is a schematic diagram of the implementation process of the object-oriented object-based non-supervised multi-temporal multispectral high-resolution remote sensing image change detection method based on CRF of the present invention;

FIG. 2(a) is a 3 rd band schematic of a high resolution IKONOS image taken in the Saudi Arabia Mina region of month 1, 2007 as used in the present invention;

FIG. 2(b) is a schematic 3 band diagram of a Mina region high resolution IKONOS image of Saudi Arabia, month 12, 2007 as used in the present invention;

fig. 2(c) is a change detection reference (Ground Truth) image;

FIG. 3(a) an EM-MRF algorithm detection result image;

FIG. 3(b) CRF algorithm detection result image;

FIG. 3(c) an image of the detection result of the OBCT algorithm;

FIG. 3(d) is an image of the results of the algorithm of the present invention.

Detailed Description

In order to make the technical means, the creation characteristics, the achievement purposes and the effects of the invention easy to understand, the invention is further described with the specific embodiments.

Referring to fig. 1, the implementation steps of the invention are as follows:

step 1: inputting two high-resolution optical remote sensing images in the same area and different time phases, and respectively recording the images as: x₁

And X₂。

Step 2: using remote sensing software ENVI to X₁And X₂Performing image registration, which comprises two steps of coarse correction and fine correction:

for geometric coarse correction, the method is realized by using related functions in ENVI4.8 software, and the specific operation steps are as follows: (1) displaying a reference image and an image to be corrected; (2) collecting ground control points GCPs; the GCPs should be uniformly distributed in the whole image, and the number of GCPs is at least more than or equal to 25. (3) Calculating an error; (4) selecting a polynomial model; (5) and (5) resampling and outputting by using bilinear interpolation.

In the bilinear difference method, when the unknown function f is found to have a value of (x, y) at the coordinate point P, it is assumed that the known function f is Q₁₁＝(x₁,y₁),Q₁₂＝(x₁,y₂),Q₂₁＝(x₂,y₁) And Q₂₂＝(x₂,y₂) Values of four coordinate points. If a coordinate system is chosen such that the coordinates of these four points are (0,0), (0,1), (1,0), and (1,1), respectively, then the bilinear interpolation formula can be expressed as:

f(x,y)≈f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy

(1)

and for geometric fine correction, performing geometric fine correction on the multispectral remote sensing image data subjected to geometric coarse correction by using an automatic matching and triangulation method.

The triangulation method comprises the steps of constructing a Delaunay triangulation network by adopting a point-by-point insertion method, determining affine transformation model parameters inside each triangle by utilizing the row and column numbers of three vertexes of each triangle and the geographic coordinates of the same-name points of the corresponding reference image, and correcting the image to be corrected to obtain the corrected remote sensing image.

And step 3: detection of X by Multivariate Change Detection (MAD) method₁And X₂And carrying out radiation normalization correction. The method comprises the steps of firstly finding out a linear combination of brightness values of each wave band of two-stage images through typical correlation analysis to obtain a difference image with enhanced change information, and obtaining a threshold value through the normalized MAD variable square sumDetermining changed and unchanged areas, and then completing relative radiation correction through a mapping equation of two-time phase pixel pairs corresponding to the unchanged areas.

And 4, step 4: and overlapping the preprocessed multi-temporal high-resolution images, and then segmenting based on an SLIC algorithm. The method comprises the following concrete steps:

4.1) initializing seed points (cluster centers): and uniformly distributing the seed points in the image according to the set number of the super pixels. Assuming that the picture has N pixel points in total and is pre-divided into R superpixels with the same size, the size of each superpixel is N/R, and the distance (step length) between adjacent seed points is approximately equal to

4.2) reselect the seed point within the nxn neighborhood of the seed point (typically taking n-3). The specific method comprises the following steps: and calculating gradient values of all pixel points in the neighborhood, and moving the seed point to the place with the minimum gradient in the neighborhood.

4.3) assigning a class label to each pixel point in the neighborhood around each seed point, wherein the search range is limited to 2S multiplied by 2S.

4.4) for each searched pixel point, respectively calculating the color distance and the space distance between the pixel point and the seed point. And taking the seed point corresponding to the minimum value as the clustering center of the pixel point.

4.5) iterate over steps 4.1-4.4 until the error converges.

4.6) solving the problems of multi-connectivity, over-small super-pixel size, cutting of a single super-pixel into a plurality of discontinuous super-pixels and the like by enhancing connectivity.

And 5: calculating the difference image of the multi-temporal image, and recording the image composed of the difference images of each wave band as X_MThen, clustering is carried out by utilizing a K-Means algorithm, and a CRF unitary energy item is determined according to a clustering result. Setting the clustering number K to 2, and implementing the following steps:

5.1) randomly selecting K clustering centers;

5.3) recalculating the centers of the classes that have been obtained.

5.4) iterating 5.2) to 5.3) until the new center is equal to the original center or less than a specified threshold value, and finishing the algorithm.

5.5) after the K-Means algorithm is converged, counting the number of objects which belong to different classes in each object, and recording the number of objects j which belong to non-variable classes as N_j1The number of the variants is N_j2Then the unary energy term for object j is:

step 6: constructing a 4-neighborhood system of CRF according to the segmentation image of the SLIC, and constructing the binary term energy of the CRF by combining the difference image in the step 5 as follows:

in the formula, i, j represents two adjacent nodes or objects of the image (in the invention, each object is regarded as one node in the graph model); c. C_i,c_jClass labels (two classes: changed and unchanged) at two nodes i, j; l (i, j) represents the length of the common boundary of the two objects; [ c ] is_i≠c_j]Denotes the 0-1 indicator function, F_M(i)、F_M(j) Respectively representing vectors consisting of the mean of the difference image values and the mean of the coordinates inside the objects i and j, F_M(i) The construction of (a) is as follows:

F_M(i)＝[Mean(Dif_i)Mean(A_ix)Mean(A_iy)] (4)

wherein Mean (Dif)_i) Representing the Mean of the difference image within the object i, Dif representing the magnitude of the difference image, Mean (A)_ix) Mean, Mean (A) representing the abscissa of all pixels within the object i_iy) Representing the mean of the ordinates of all the pixels within the object i.

The final CRF energy term constructed was as follows:

in the formula of U_j、U_i,jA unitary energy term representing an object j or a node j (in the present invention, each object is regarded as a node in the graph model), a binary energy term between adjacent objects j and i, S represents a set of nodes, E represents a set of adjacent nodes, and a parameter β is a regularization parameter for controlling a specific gravity between the unitary energy term and the binary energy term.

Step 10: and optimizing the CRF energy item U according to an optimization algorithm cycle Belief Propagation (LBP) algorithm to obtain a final change detection result.

The experimental data of the invention are multi-temporal IKONOS high-resolution image data in the Mina region of Saudi Arabia, the image size is 700 x 950, and three wave bands of B1, B2 and B3 are used. To verify the effectiveness of the present invention, the change detection method of the present invention was compared to the following change detection methods:

(1) CVA-based EM-MRF method (EM-MRF) [ detection methods mentioned in the article "Automatic analysis of differential image for unsupervised change detection" (IEEE Transactions on Geoscience and Remote Sensing,2000,38(3):1171-1182 ], by Bruzzone L. et al, Italy ].

(2) The CRF detection (CRF) method proposed by Cao et al [ Guo Cao, Xuesong Li & Liuun Zhou. Un. upstream changed detection in high spatial resolution removal Sensing fields on a conditional random field model. European Journal of removal Sensing,2016,49:225-

(3) Object-oriented chi-square transform-based detection method (OBCT) proposed by Descl e et al [ B.Descl e, P.Bogart, and P.Defourny, "Forest change detection by static object-based method," Remote Sens. environ,2006,102 (1-2), 1-11 ]

(4) The method of the invention.

The detection performance is measured by four indexes of error detection number FP, missing detection number FN, total error number OE and Kappa coefficient. The closer FP, FN and OE are to 0 and the closer Kappa coefficient is to 1, indicating the better performance of the change detection method. The results are shown in Table 1.

TABLE 1 comparison of multi-temporal IKONOS image change detection results in Mina area

Method	FP	FN	OE	k
					EM-MRF	13373	1367	14740	0.821
CRF	650	3913	4563	0.934
					OBCT	8813	780	9593	0.878
The method of the invention	1624	2139	3763	0.947
					Ideal for	0	0	0	1

As can be seen from Table 1, the Kappa coefficient of the detection method provided by the invention is the largest and is closer to 1 than the other three detection algorithms. In addition, the total error number OE of the present invention is the smallest in the comparison algorithm, closer to 0. In conclusion, the performance of the change detection algorithm of the invention is superior to that of the other three detection methods, which shows that the change detection method provided by the invention is effective.

Fig. 2(a) is a previous-time-phase multispectral IKONOS image of the Mina region, fig. 2(b) is a next-time-phase multispectral IKONOS image of the Mina region, and fig. 2(c) is a reference image of change detection. Fig. 3(a) is a change detection result of the EM-MRF algorithm, fig. 3(b) is a change detection result of the CRF algorithm, fig. 3(c) is a change detection result of the OBCT algorithm, and fig. 3(d) is a change detection result of the method of the present invention. From the comparison between the reference diagram of fig. 2(c) and fig. 3(a) - (d), the detection effect of the algorithm of the present invention is the best in terms of visual effect.

The foregoing shows and describes the general principles and broad features of the present invention and advantages thereof. It will be understood by those skilled in the art that the present invention is not limited to the embodiments described above, which are described in the specification and illustrated only to illustrate the principle of the present invention, but that various changes and modifications may be made therein without departing from the spirit and scope of the present invention, which fall within the scope of the invention as claimed. The scope of the invention is defined by the appended claims and equivalents thereof.

Claims

1. An object-oriented multispectral high-resolution remote sensing image change detection method is characterized by comprising the following steps:

and 5: calculating a difference image of the multi-temporal images, and recording the difference image as X_MThen, clustering by using a K-Means algorithm, and determining a CRF unitary energy item according to a clustering result;

step 6: constructing a neighborhood system of the CRF according to the segmentation image of the SLIC, and constructing binary item energy of the CRF according to the difference image and the space coordinate of the object in the step 5;

and 7: and optimizing the CRF energy item U according to a circulation reliability propagation (LBP) optimization algorithm to obtain a final change detection result.

2. The method according to claim 1, wherein the step 2 of coarse correction comprises the following specific steps:

2.1) displaying the reference image and the image to be corrected;

2.3) calculating the root mean square error;

2.4) selecting a polynomial correction model;

2.5) resampling and outputting by using bilinear interpolation.

3. The method according to claim 2, wherein the unknown function f is evaluated at a coordinate point P (x, y), x and y respectively representing an abscissa and an ordinate; based on a predetermined known function f at Q₁₁＝(x₁,y₁),Q₁₂＝(x₁,y₂),Q₂₁＝(x₂,y₁) And Q₂₂＝(x₂,y₂) Values of four coordinate points, x₁，x₂Is the abscissa, y₁，y₂Is ordinate, if a coordinate system is selected such that the coordinates of these four points are (0,0), (0,1), (1,0), and (1,1), respectively, then the bilinear interpolation formula is expressed as:

f(x,y)≈f(0,0)(1-x)(1-y)+f(1,0)x(1-y)+f(0,1)(1-x)y+f(1,1)xy (1) 。

4. the method according to claim 1, wherein in step 2, the fine correction is performed on the multispectral remote sensing image data to be geometrically coarsely corrected by using an automatic matching and triangulation method.

5. The method according to claim 4, wherein the triangulation method is a Delaunay triangulation method constructed by a point-by-point interpolation method, for each triangle, affine transformation model parameters inside the triangle are determined by using the row and column numbers of three vertexes of the triangle and the geographic coordinates of the same-name points of the corresponding reference image, and the image to be corrected is corrected to obtain the corrected remote sensing image.

6. The object-oriented multispectral high-resolution remote sensing image change detection method according to claim 1, wherein in step 3, the specific method for performing normalization correction on the MAD relative to radiation in the multivariate change detection is as follows:

7. The object-oriented multispectral high-resolution remote sensing image change detection method according to claim 1, wherein the step 4 is implemented by the following steps:

4.1) initializing seed points, namely cluster centers: uniformly distributing seed points in the image according to the number of the preset super pixels; based on the total N pixel points of the preset image, the image is pre-divided into R super-pixels with the same size, the size of each super-pixel is N/R, and the distance between adjacent seed points, namely the step length, is

4.5) iterating step 4.1) -step 4.4) until the error converges.

8. The object-oriented multispectral high-resolution remote sensing image change detection method according to claim 7, wherein in step 4.2), the specific selection method of the seed points is as follows:

9. The object-oriented multispectral high-resolution remote sensing image change detection method according to claim 1, wherein in step 5, based on a preset cluster number K being 2, the CRF unitary energy item is specifically determined by the following method:

5.1) randomly selecting K clustering centers;

5.3) recalculating centers of the obtained classes;

10. the object-oriented multispectral high-resolution remote sensing image change detection method according to claim 1, wherein in step 6, the binary term energy for constructing the CRF is as follows:

in the formula, i, j represents two adjacent nodes of the image, namely adjacent objects; c. C_i,c_jRepresenting class labels at two nodes i, j; two categories are changed and unchanged; l (i, j) represents the length of the common boundary of the two objects; [ c ] is_i≠c_j]Represents a 0-1 indicator function; f_M(i)、F_M(j) Respectively representing vectors consisting of the mean of the difference image values and the mean of the coordinates inside the objects i and j, wherein F_M(i) The construction of (a) is as follows:

F_M(i)＝[Mean(Dif_i) Mean(A_ix) Mean(A_iy)] (4)

the final CRF energy term constructed was as follows:

wherein U represents the total energy term for CRF, U_j、U_i,jThe parameter beta is a regularization parameter used for controlling the proportion between the unary energy item and the binary energy item.