CN110276746B - Robust remote sensing image change detection method - Google Patents

Abstract

The invention discloses a robustness remote sensing image change detection method, which comprises the following steps: acquiring image information of two time phases of a target area; preprocessing the remote sensing images of the two time phases to obtain image characteristics; selecting and screening image features according to a Relief algorithm; performing feature dimension reduction by using a PCA algorithm to obtain a final preferred feature change vector; and (4) classifying by using an NR-Kmeans method to obtain a change detection result. The invention improves the traditional kmeans clustering algorithm, and has higher detection precision and high accuracy.

Description

A Robust Remote Sensing Image Change Detection Method

technical field

The invention belongs to the technical field of image changes, and in particular relates to a robust remote sensing image change detection method.

Background technique

Remote sensing, that is, long-distance perception. The term remote sensing first appeared in a military scientific research program in the United States in 1960, and was formally adopted in the first academic seminar on environmental remote sensing held by the University of Michigan and other units in 1962. Nowadays, remote sensing technology is used in more and more fields, such as: dynamic change monitoring of forest or vegetation, post-disaster analysis and assessment of natural disasters, analysis of land use changes, monitoring of farmland, real-time monitoring of urban changes, Analyzing the growth status of crops and conducting dynamic surveillance of military strategic targets, such as airports and roads, has greatly promoted economic and social development.

Remote sensing image change detection belongs to the field of remote sensing image processing, and is used to analyze and process remote sensing images of the same location in different periods to obtain change information. Scholars at home and abroad have carried out a lot of research on the problem of remote sensing image change detection, and proposed various change detection methods. Among them, the feature-based change detection method is to use the spectral features, texture features, and spatial structure features in the image to detect changes in remote sensing images. The feature-based change detection method considers the rich feature information of the image, and has good noise robustness and anti-light and radiation interference ability.

Since remote sensing images are affected by illumination, atmospheric environment, etc., the collected image features are also affected by them. Not only that, due to noise interference, the traditional clustering algorithm is affected by sparse data and the classification results are not accurate enough. Selecting representative features and accurately classifying them becomes the key to change detection technology in remote sensing images.

SUMMARY OF THE INVENTION

Purpose of the invention: In view of the above problems, the present invention proposes a robust remote sensing image change detection method. This method is a feature-based remote sensing image change detection method. Image features are selected and screened according to the Relief algorithm, then the PCA algorithm is used to reduce the dimension, and finally the NR-Kmeans algorithm is used to eliminate the influence of noise interference points on the detection results. The final output result. The invention eliminates noise interference points, and the detection is more accurate.

Technical solution: In order to achieve the purpose of the present invention, the technical solution adopted in the present invention is: a robust remote sensing image change detection method, comprising the following steps:

( ₁ ) collect and obtain remote sensing images of target areas T1 and _T2 , which meet the resolution requirements and contain clear texture feature information;

(2) Using ENVI software to realize image data preprocessing to eliminate shooting errors; using the eCogntion software platform to perform multi-scale segmentation on the preprocessed images at T ₁ and T ₂ , and divide the two images into N image blocks respectively, That is, N objects are obtained respectively; the spectral features and texture features of the images at times T ₁ and T ₂ are obtained respectively, and feature vectors are constructed to obtain a change vector set M _d , where M _d =[m ₁ ,m ₂ ,...m _N ], wherein, m _i represents the change vector corresponding to the ith object; the elements in M _d are arranged according to the arrangement of N objects in the image at time T ₁ or T ₂ to form a difference map;

(3) According to the Relief feature selection algorithm, obtain the weight of the feature of each object in the feature subspace M _l of the change vector set M _d , select e change features with the largest weight, and the value of e is less than the total number of features contained in the object, and then The e change features are selected from the N objects respectively, and the matrix X is obtained in combination with the change vector set M _d ;

(4) Use PCA algorithm to reduce the feature dimension of matrix X, and obtain the changed feature matrix Y after transformation, Y=[y ₁ , y ₂ ,...,y _N ];

(5) Use the NR-Kmeans algorithm to classify the elements _yi of the variable feature matrix Y, remove the sparse data disturbed by noise from the variable feature matrix Y, and obtain a dense point set Y';

(6) Divide the dense point set Y′ into the unchanged class and the changed class. The gray value of the object in the difference map corresponding to the unchanged class element is converted to 255, and the difference value corresponding to the element in the changed class The gray value of the object in the figure is converted to 0, and the change detection result map is output according to the gray value of the object, and the area of the image change is marked in the figure.

Further, in step (2), ENVI software is used to realize image data preprocessing to eliminate shooting errors; eCogntion software platform is used to perform multi _- scale segmentation on the preprocessed images at T1 and T2, and the _two images are divided into N image blocks, that is, N objects are obtained respectively; the spectral features and texture features of the images at time T ₁ and T ₂ are obtained respectively, and feature vectors are constructed to obtain a set of change vectors. The steps are as follows:

( _2.1 ) Using ENVI software to perform geometric correction, image registration, and atmospheric correction operations _on the images at T1 and T2, respectively, to realize image data preprocessing;

(2.2) Use the eCogntion software platform to perform multi-scale segmentation on the preprocessed images at T ₁ and T ₂ , and divide the two images into N image blocks respectively, that is, to obtain N objects respectively;

(2.3) Obtain the Mean-Std feature of the preprocessed image as the spectral feature of the image; use the gray level co-occurrence matrix to extract the texture feature of the preprocessed image;

(2.4) According to the spectral features and texture features extracted from the two preprocessed images, a feature vector is constructed, and a change vector set M _d is obtained.

Further, the Mean-Std feature of the image is obtained as described in step (2.3) as the spectral feature of the image; the method is as follows:

For each object of the image, the mean and standard deviation Std in the spectral features are obtained according to the following formulas:

In the formula, Mean is the mean of the grayscale of the pixels in the object, std is the standard deviation of the grayscale of the pixels in the object, A represents the number of pixels in the object, and c _i represents the grayscale of the ith pixel in the object.

Further, in step (2.3), the gray level co-occurrence matrix is used to extract the texture features of the image; the method is as follows:

(2.3.1) Any pixel point (τ ₁ , τ ₂ ) and another pixel point (τ ₁ +Δ ₁ , τ ₂ +Δ ₂ ) in the image form a point pair, where Δ ₁ , Δ ₂ are integers; Assuming that the gray value of the pixel point (τ ₁ , τ ₂ ) is f1, the gray value of (τ ₁ +Δ ₁ , τ ₂ +Δ ₂ ) is f2, then the gray value of the pair of points is (f1, f2 ), set the maximum gray level of the image to be L, then the combination of f1 and f2 has a total of L*L groups;

(2.3.2) Count the number of occurrences of each group of different (f1, f2) values in the image, and then arrange them into a matrix, where the matrix element values located at the matrix position (L, L) are the two grayscale values. are the number of occurrences of point pairs of L;

(2.3.3) According to the total number of occurrences of (f1, f2), the number of occurrences of each group of different (f1, f2) values is normalized to the probability of occurrence g(f1, f2), and the normalized The square matrix is called the gray level co-occurrence matrix;

(2.3.4) The grayscale co-occurrence matrix provides 14 kinds of statistics as texture features, the present invention selects 6 kinds of statistics as texture features, and the 6 kinds of statistics are contrast, entropy, energy, equality, dissimilarity, Correlation; contrast reflects the clarity of the image and the depth of texture grooves; entropy expresses the degree of non-uniformity or complexity of the texture in the image; energy reflects the uniformity of the grayscale distribution of the image and the thickness of the texture; equality reflects the image The homogeneity of the texture measures the local variation of the image texture; the correlation reflects the consistency of the image texture; these 6 kinds of statistics represent the characteristics of the image from various aspects, so these 6 kinds of statistics are selected as texture features;

(2.3.5) Assuming that the value of the matrix element at the position (i, j) of the gray level co-occurrence matrix is g(i, j), and L is the maximum gray level of the image, the texture feature is expressed as follows:

In the formula, Con represents contrast, Ent represents entropy, Ene represents energy, Hom represents equality, Dis represents dissimilarity, and Cor represents correlation;

Further, according to the spectral features and texture features extracted from the two preprocessed images described in step (2.4), a change vector set M _d is obtained; the method is as follows:

(2.4.1) According to the spectral features and texture features of the images at time T ₁ and T ₂ obtained in step (2.3), construct feature vectors respectively, and calculate the difference between the corresponding feature vectors of the two images; the two images correspond to The difference M _d (i, j) of the j-th eigenvector of the i-th object is expressed as follows:

M _d (i, j) = M _{d, 1} (i, j) - M _{d, 2} (i, j)

In the formula, M _{d, 1} (i, j) is the j-th feature vector of the i-th object in the image at time T ₁ , and M _{d, 2} (i, j) is the j-th object of the i-th object in the image at time T ₂ feature vector; i≤N, N is the total number of objects in the image; j≤S _feature , S _feature is the total number of features contained in the object; S _feature =S _band ×8, S _band is the total number of bands of the object;

According to the difference value M _d (i, j), the change vector m _i corresponding to the ith object is obtained, which is expressed as follows:

m _i = (M _d (i, 1), M _d (i, 2), ..., M _d (i, S _feature ))

(2.4.2) Repeat step (2.4.1), calculate the difference between the feature vectors of each object corresponding to the two images, obtain the change vector corresponding to each object, and form the change vector set M _d , which is expressed as follows :

M _d =[m ₁ , m ₂ , ... m _N ]

Arrange the elements in M _d according to the arrangement of N objects in the image at time T1 and _T2 to form _a difference map.

Further, according to the Relief feature selection algorithm described in step (3), obtain the weight of the feature of each object in the feature subspace M _i of the change vector set M _d , select e the change features with the largest weight, and the value of e is less than the object contains The total number of features, and then select the e change features from the N objects, and combine the change vector set M _d to obtain the matrix X; the method is as follows:

(3.1) The change vector set M _d is also called the feature space, and the subspace M _l =[m ₁ ,m ₂ ,m ₃ ] of the feature space M _d =[m ₁ ,m ₂ ,m ₃ ,... m _N ] is extracted, ... m _n ] as the training set, where n<N; the feature subspace M _l is divided into two number sets M _l1 and M _l2 using artificial interpretation, where M _l1 contains l1 elements, M _L2 contains l2 elements, l1+l2=n, M _l1 ∈ M _l , M _l2 ∈ M _l ;

(3.2) Calculate the Euclidean distance of the change vector between all elements in M _l to form a distance matrix, that is

表示M_l中元素两两之间变化矢量的欧式距离；where d represents the distance matrix,

represents the Euclidean distance of the change vector between the elements in M _l ;

(3.3) In the feature subspace M _l , for the change vector m _i of the object i, i=1, 2, 3, ..., n, if m _i ∈ M _l1 , select the one with the smallest Euclidean distance in M _l1 s vectors, and select the s vectors with the largest Euclidean distance in M _l2 ; if m _i ∈ M _l2 , select the s vectors with the smallest Euclidean distance in M _l2 , and select the largest Euclidean distance in M _l1 s vectors; that is, the nearest neighbor matrix Q _i = [q ₁ , q ₂ ,..., q _s ,..., q _2s ] for _mi , where [q ₁ , q ₂ ,..., q _s ] are the s vectors with the smallest Euclidean distance, [q _s+1 , q _s+2 , ..., q _2s ] are the s vectors with the largest Euclidean distance, s=min(l1, l2);

(3.4) Calculate the weight corresponding to each feature of object i according to the related statistics of the wrongly guessed neighbors and the related statistics of the correctly guessed neighbors, namely

是特征f_j对应的权重，diffr_j是特征f_j的猜错近邻相关统计量，diffw_j是特征f_j的猜中近邻相关统计量；diffr_j和diffw_j的计算表达式如下：In the formula,

is the weight corresponding to the feature f _j , diffr _j is the related statistic of the wrong-guessed neighbor of the feature f _j , diffw _j is the related statistic of the correct-guessed neighbor of the feature f _j ; the calculation expressions of diffr _j and diffw _j are as follows:

In the formula, s=min(l1, l2), i represents the object i, j represents the jth feature, when a=1, 2, 3,..., s, Q _i (a, j) represents [q _{1 , q 2 , .} ..., q _2s ]; when b = 1, 2, 3, ..., s, Q _i (b, j) represents [q ₁ , q ₂ , ..., q _s ], when b = When s+1, s+2, ..., 2s, Q _i (b, j) represents [q _s+1 , q _s+2 , ..., q _2s ];

(3.5) According to the weight of the feature of each object in the feature subspace M _l calculated in step (3.4), select e change features with the largest weight, and then select the e change features from the N objects respectively, combine with From the original change vector set M _d , the matrix X is obtained, which is expressed as follows:

In the formula, X _i =[x _i1 x _i2 ... x _ie ], i=1, 2, ..., N; the matrix element x _ij represents the j-th feature in the e changing features of the object i, e The value is less than the total number of features contained in the object.

Further, described in step (4), the PCA algorithm is used to perform feature dimension reduction on the matrix X to obtain the change feature matrix Y after the transformation; the method is as follows:

(4.1) According to the matrix X obtained in step (3), calculate the mean value μ _i of each row of the matrix, i=1, 2,...,N:

In the formula, j=1, 2,...,e, E(X _i ) represents mathematical expectation;

(4.2) Calculate the covariance matrix C of the matrix X according to the mean μ _i , namely:

In the formula, i=1,2,...,N,r=1,2,...,N,j=1,2,...,e; it can be seen from the above formula that the covariance matrix C is a real Symmetric matrix;

(4.3) Calculate the eigenvalues λ ₁ , λ ₂ , ..., λ _N of the covariance matrix C, and the corresponding eigenvectors are v ₁ , v ₂ , ..., v _N , that is, λ _i v _i =Cv _i , i=1,2,...,N;

(4.4) Sort the eigenvalues λ ₁ , λ ₂ , ..., λ _N in descending order to obtain λ' ₁ , λ' ₂ , ..., λ' _N , and the corresponding eigenvectors are v' ₁ , v' ₂ ,...,v′ _N , take the sorted eigenvalues to form a diagonal matrix P;

(4.5) Calculate the projection of the matrix X on the new eigenvectors v′ ₁ , v′ ₂ , ..., v′ _N , and obtain the changed feature matrix Y after dimension reduction, namely:

Y=X*P=[y ₁ , y ₂ , ..., y _N ]

Among them, y ₁ , y ₂ , ..., y _N represent the changing characteristics of the object.

Further, in the step (5), the NR-Kmeans algorithm is used to classify the element _yi of the variable feature matrix Y, and the sparse data disturbed by noise is removed from the variable feature matrix Y to obtain a dense point set Y'; the steps are as follows:

的元素点

Q_b(y_α)表示包含元素点y_α在内的y_α的最近邻元素点集合；其中，

表示元素点y_α、

的距离，b为y_α在邻域范围内最近邻元素点个数，f＝1，2，...，b；(5.1) In the matrix Y=[y ₁ , y ₂ , ..., y _N ], for each element point y _α , y _α ∈ Y, given the neighborhood radius δ, get the

element point of

Q _b (y _α ) represents the nearest neighbor element point set of y _α including element point y _α ; where,

represents the element point y _α ,

distance, b is the number of nearest neighbor element points of y _α in the neighborhood, f=1, 2, ..., b;

(5.2) Calculate the density function value of the element point y _α ; the density function value of the element point y _α represents the sum of the functions of all the nearest neighbor element points in the range of the neighborhood radius δ; the calculation formula of the density function value of y _α is as follows :

为元素点y_α、

的距离，f＝1，2，...，b；In the formula,

is the element point y _α ,

distance, f = 1, 2, ..., b;

(5.3) Repeat steps (5.1) to (5.2) until the calculation of the density function values of all element points in the matrix Y is completed;

(5.4) When the density function value of y _β is greater than or equal to the average density value of all element points in the matrix Y, the element point y _β is regarded as available data and put into the dense point set Y′, that is, y _β conforms to:

In the formula, DF(y _β ) represents the density function value of y _β , and DF(y _α ) is the density function value of the data point y _α ;

(5.5) Obtain all available data in the matrix Y according to step (5.4), the available data constitute a dense point set Y', and Y' is an element point set after excluding sparse data interfered by noise.

Further, according to step (6), the dense point set Y' is divided into unchanged classes and changed classes, and the gray value of the objects in the difference map corresponding to the unchanged class elements is converted to 255, and the changed classes The gray value of the object in the difference map corresponding to the element is converted to 0, and the change detection result map is output according to the gray value of the object; the method is as follows:

(6.1) From the dense point set Y′=[y′ ₁ , y′ ₂ , ..., y′ _η ], η is the number of element points in Y′, η≤N, y′ _η means that the noise is eliminated The characteristics of the object behind the disturbed data; select the element point with the largest density value as the first initial cluster center C ₁ , and calculate the difference between the eigenvalues of other element points in Y′ and the eigenvalues of the cluster center C ₁ , select the element point with the largest difference from the _{eigenvalue of} C ₁ , that is, the point with the largest absolute value of the difference as the _second initial cluster center C _{2 ;} Class centers C ₁ and C ₂ correspond to classes D ₁ and D ₂ ;

(6.2) Calculate the difference between the element point y′ _θ in the dense point set Y′ and the eigenvalues of the two cluster centers, namely

d _θ1 = |y′ _θ -c ₁ |

d _θ2 = |y′ _θ -c ₂ |

d _θ1 represents the absolute value of the feature difference between the element point y' _θ and the cluster center C ₁ , d _θ2 represents the absolute value of the feature difference between the element point y' _θ and the cluster center C ₂ , θ=1, 2, ..., η; divide the elements in the dense point set Y' into the class corresponding to the cluster center with the smaller absolute difference value;

(6.3) Calculate the average value c′ ₁ , c′ ₂ of the corresponding element point eigenvalues in the classes D ₁ , D ₂ , namely

In the formula, c′ ₁ , c′ ₂ are the average values of the corresponding element point eigenvalues of class D ₁ and D ₂ respectively; |D ₁ | represents the number of element points in class D ₁ , and |D ₂ | represents class D The number of element points in ₂ , y′ corresponds to the eigenvalues of the element points in the class; if c′ ₁ ≠c ₁ , update the eigenvalue c ₁ of the cluster center C ₁ of class D ₁ to c′ ₁ ; if c′ 1 ≠c 1 , ₂ ≠c ₂ , update the eigenvalue c ₂ of the cluster center C ₂ of class D ₂ to c′ ₂ ;

(6.4) Repeat steps (6.2) to (6.3) until the eigenvalues of the two cluster centers no longer change;

(6.5) Divide the two classes D ₁ , D ₂ into the unchanged class and the changed class, wherein the class with a smaller average value is the unchanged class, and the class with a larger average value is the changed class; the unchanged class The gray value of the object in the difference map corresponding to the element is converted to 255, and the gray value of the object in the difference map corresponding to the element in the changed class is converted to 0;

(6.6) Output the detection result graph according to the gray value of the object, and the area where the image changes is marked in the graph.

Beneficial effects: compared with the prior art, the technical solution of the present invention has the following beneficial technical effects:

(1) In order to reduce the interference of noise on the data, the density distribution function is used to eliminate the density sparse points and some noise interference points.

(2) Combined with the principle of the farthest distance, the initial cluster center is selected in the high-density data set, which avoids the disadvantage of randomly selecting the initial cluster center in the original Kmeans algorithm.

Description of drawings

Fig. 1 is the flow chart of the present invention;

Figure 2 is a flowchart of the NR-Kmeans algorithm.

Detailed ways

The technical solutions of the present invention will be further described below with reference to the accompanying drawings and embodiments.

A robust remote sensing image change detection method according to the present invention, as shown in Figure 1, includes the following steps:

( ₁ ) Collect and obtain the remote sensing images of the target areas T1 and _T2 , which meet the resolution requirements and contain clear texture feature information; select the image of a certain area in Hefei City, including the coverage of farmland crops;

(2) Using ENVI software to realize image data preprocessing to eliminate shooting errors; using the eCogntion software platform to perform multi-scale segmentation on the preprocessed images at T ₁ and T ₂ , and divide the two images into N image blocks respectively, That is, N objects are obtained respectively; the spectral features and texture features of the images at times T ₁ and T ₂ are obtained respectively, and feature vectors are constructed to obtain a change vector set M _d , where M _d =[m ₁ ,m ₂ ,...m _N ], wherein, m _i represents the change vector corresponding to the ith object; the elements in M _d are arranged according to the arrangement of N objects in the image at time T ₁ or T ₂ to form a difference map;

(3) According to the Relief feature selection algorithm, obtain the weight of the feature of each object in the feature subspace M _l of the change vector set M _d , select e change features with the largest weight, and the value of e is less than the total number of features contained in the object, and then The e change features are selected from the N objects respectively, and the matrix X is obtained in combination with the change vector set M _d ;

(4) Use the PCA algorithm to reduce the feature dimension of the matrix X, and obtain the changed feature matrix Y after the transformation, Y=[y ₁ , y ₂ ,...,y _N ];

(5) Use the NR-Kmeans algorithm to classify the elements _yi of the variable feature matrix Y, remove the sparse data disturbed by noise from the variable feature matrix Y, and obtain a dense point set Y';

(6) Divide the dense point set Y' into the unchanged class and the changed class. The gray value of the object in the difference map corresponding to the unchanged class element is converted to 255, and the difference value corresponding to the element in the changed class The gray value of the object in the figure is converted to 0, and the change detection result map is output according to the gray value of the object, and the area of the image change is marked in the figure.

In step (2), ENVI software is used to realize image data preprocessing to eliminate shooting errors; eCogntion software platform is used to perform multi _- scale segmentation on the preprocessed images at time T1 and T2, and the _two images are divided into N pieces respectively. Image blocks, that is, N objects are obtained respectively; the spectral features and texture features of the images at times T ₁ and T ₂ are obtained respectively, and feature vectors are constructed to obtain a set of change vectors. The steps are as follows:

(2.1) Using ENVI software to perform geometric correction, image registration, and atmospheric correction operations on the images at T ₁ and T ₂ , respectively, to realize image data preprocessing; the size of both images is 512*512;

(2.2) Use the eCogntion software platform to perform multi-scale segmentation on the preprocessed images at T ₁ and T ₂ , and divide the two images into N image blocks respectively, that is, to obtain N objects respectively;

(2.3) Obtain the Mean-Std feature of the preprocessed image as the spectral feature of the image; use the gray level co-occurrence matrix to extract the texture feature of the preprocessed image;

(2.4) According to the spectral features and texture features extracted from the two preprocessed images, a feature vector is constructed, and a change vector set M _d is obtained.

The Mean-Std feature of the image is obtained as described in step (2.3) as the spectral feature of the image; the method is as follows:

For each object of the image, the mean and standard deviation Std in the spectral features are obtained according to the following formulas:

In the formula, Mean is the mean of the grayscale of the pixels in the object, std is the standard deviation of the grayscale of the pixels in the object, A represents the number of pixels in the object, and c _i represents the grayscale of the ith pixel in the object.

In step (2.3), the texture feature of the image is extracted by using the grayscale co-occurrence matrix; the method is as follows:

(2.3.1) Any pixel point (τ ₁ , τ ₂ ) and another pixel point (τ ₁ +Δ ₁ , τ ₂ +Δ ₂ ) in the image form a point pair, where Δ ₁ , Δ ₂ are integers; Assuming that the gray value of the pixel point (τ ₁ , τ ₂ ) is f1, the gray value of (τ ₁ +Δ ₁ , τ ₂ +Δ ₂ ) is f2, then the gray value of the pair of points is (f1, f2 ), set the maximum gray level of the image to be L, then the combination of f1 and f2 has a total of L*L groups; in this embodiment, if L=16, then the combination of f1 and f2 has a total of 16*16 groups;

(2.3.2) Count the number of occurrences of each group of different (f1, f2) values in the image, and then arrange them into a matrix, where the matrix element values located at the matrix position (L, L) are the two grayscale values. are the number of occurrences of point pairs of L;

(2.3.3) According to the total number of occurrences of (f1, f2), the number of occurrences of each group of different (f1, f2) values is normalized to the probability of occurrence g(f1, f2), and the normalized The square matrix is called the gray level co-occurrence matrix;

(2.3.4) The grayscale co-occurrence matrix provides 14 kinds of statistics as texture features, the present invention selects 6 kinds of statistics as texture features, and the 6 kinds of statistics are contrast, entropy, energy, equality, dissimilarity, Correlation; contrast reflects the clarity of the image and the depth of texture grooves; entropy expresses the degree of non-uniformity or complexity of the texture in the image; energy reflects the uniformity of the grayscale distribution of the image and the thickness of the texture; equality reflects the image The homogeneity of the texture measures the local variation of the image texture; the correlation reflects the consistency of the image texture; these 6 kinds of statistics represent the characteristics of the image from various aspects, so these 6 kinds of statistics are selected as texture features;

(2.3.5) Assuming that the value of the matrix element at the position (i, j) of the gray level co-occurrence matrix is g(i, j), and L is the maximum gray level of the image, the texture feature is expressed as follows:

In the formula, Con represents contrast, Ent represents entropy, Ene represents energy, Hom represents equality, Dis represents dissimilarity, and Cor represents correlation;

Described in step (2.4), according to the spectral features and texture features extracted from the two preprocessed images, obtain the change vector set M _d ; the method is as follows:

(2.4.1) According to the spectral features and texture features of the images at time T ₁ and T ₂ obtained in step (2.3), construct feature vectors respectively, and calculate the difference between the corresponding feature vectors of the two images; the two images correspond to The difference M _d (i, j) of the j-th eigenvector of the i-th object is expressed as follows:

M _d (i, j) = M _{d, 1} (i, j) - M _{d, 2} (i, j)

In the formula, M _{d, 1} (i, j) is the j-th feature vector of the i-th object in the image at time T ₁ , and M _{d, 2} (i, j) is the j-th object of the i-th object in the image at time T ₂ feature vector; i≤N, N is the total number of objects in the image; j≤S _feature , S _feature is the total number of features included in the object; S _feature =S _band × 8, S _band is the total number of bands of the object; in this embodiment, the image There are 4 bands, then S _feature = 32;

According to the difference value M _d (i, j), the change vector m _i corresponding to the ith object is obtained, which is expressed as follows:

m _i = (M _d (i, 1), M _d (i, 2), ..., M _d (i, 32))

(2.4.2) Repeat step (2.4.1), calculate the difference between the feature vectors of each object corresponding to the two images, obtain the change vector corresponding to each object, and form the change vector set M _d , which is expressed as follows :

M _d =[m ₁ , m ₂ , . . . m _N ]

Arrange the elements in M _d according to the arrangement of N objects in the image at time T1 and _T2 to form _a difference map.

According to the Relief feature selection algorithm described in step (3), the weight of the feature of each object in the feature subspace M _l of the change vector set M _d is obtained, and the e change features with the largest weight are selected, and the value of e is less than the feature contained in the object. the total number, and then select the e change features from the N objects respectively, and combine the change vector set M _d to obtain the matrix X; the method is as follows:

(3.1) The change vector set M _d is also called the feature space, extract the subspace M _l =[m ₁ ,m ₂ ,m of the feature space M _d =[m ₁ ,m ₂ ,m ₃ ,... m _N ] _{3 , ... _ _ _ _} l2 elements, l1+l2=n, M _l1 ∈ M _l , M _l2 ∈ M _l ;

(3.2) Calculate the Euclidean distance of the change vector between all elements in M _l to form a distance matrix, that is

表示M_l中元素两两之间变化矢量的欧式距离；where d represents the distance matrix,

represents the Euclidean distance of the change vector between the elements in M _l ;

(3.3) In the feature subspace M _l , for the change vector m _i of the object i, i=1, 2, 3, ..., n, if m _i ∈ M _l1 , select the one with the smallest Euclidean distance in M _l1 s vectors, and select the s vectors with the largest Euclidean distance in M _l2 ; if m _i ∈ M _l2 , select the s vectors with the smallest Euclidean distance in M _l2 , and select the largest Euclidean distance in M _l1 s vectors; that is, the nearest neighbor matrix Q _i = [q ₁ , q ₂ ,..., q _s ,..., q _2s ] for _mi , where [q ₁ , q ₂ ,..., q _s ] are the s vectors with the smallest Euclidean distance, [q _s+1 , q _s+2 , ..., q _2s ] are the s vectors with the largest Euclidean distance, s=min(l1, l2);

(3.4) Calculate the weight corresponding to each feature of object i according to the related statistics of the wrongly guessed neighbors and the related statistics of the correctly guessed neighbors, namely

是特征f_j对应的权重，diffr_j是特征f_j的猜错近邻相关统计量，diffw_j是特征f_j的猜中近邻相关统计量；diffr_j和diffw_j的计算表达式如下：In the formula,

is the weight corresponding to the feature f _j , diffr _j is the related statistic of the wrong-guessed neighbor of the feature f _j , diffw _j is the related statistic of the correct-guessed neighbor of the feature f _j ; the calculation expressions of diffr _j and diffw _j are as follows:

In the formula, s=min(l1, l2), i represents the object i, j represents the jth feature, when a=1, 2, 3,..., s, Q _i (a, j) represents [q _{1 , q 2 , .} ..., q _2s ]; when b = 1, 2, 3, ..., s, Q _i (b, j) represents [q ₁ , q ₂ , ..., q _s ], when b = When s+1, s+2, ..., 2s, Q _i (b, j) represents [q _s+1 , q _s+2 , ..., q _2s ];

(3.5) According to the weight of the feature of each object in the feature subspace M _l calculated in step (3.4), select e change features with the largest weight, and then select the e change features from the N objects respectively, and combine From the original change vector set M _d , the matrix X is obtained, which is expressed as follows:

In the formula, X _i =[x _i1 x _i2 ... x _ie ], i=1, 2, ..., N; the matrix element x _ij represents the j-th feature in the e changing features of the object i, e The value is less than the total number of features contained in the object.

The described step (4) utilizes PCA algorithm to carry out feature dimension reduction to matrix X, and obtains the change feature matrix Y after the transformation; The method is as follows:

(4.1) According to the matrix X obtained in step (3), calculate the mean value μ _i of each row of the matrix, i=1, 2,...,N:

In the formula, j=1, 2,...,e, E(X _i ) represents mathematical expectation;

(4.2) Calculate the covariance matrix C of the matrix X according to the mean μ _i , namely:

In the formula, i=1,2,...,N,r=1,2,...,N,j=1,2,...,e; it can be seen from the above formula that the covariance matrix C is a real Symmetric matrix;

(4.3) Calculate the eigenvalues λ ₁ , λ ₂ , ..., λ _N of the covariance matrix C, and the corresponding eigenvectors are v ₁ , v ₂ , ..., v _N , that is, λ _i v _i =Cv _i , i=1, 2, ..., N;

(4.4) Sort the eigenvalues λ ₁ , λ ₂ , ..., λ _N in descending order to obtain λ' ₁ , λ' ₂ , ..., λ' _N , and the corresponding eigenvectors are v' ₁ , v' ₂ ,...,v′ _N , take the sorted eigenvalues to form a diagonal matrix P;

(4.5) Calculate the projection of the matrix X on the new eigenvectors v′ ₁ , v′ ₂ , ..., v′ _N , and obtain the changed feature matrix Y after dimension reduction, namely:

Y=X*P=[y ₁ , y ₂ , ..., y _N ]

Among them, y ₁ , y ₂ , ..., y _N represent the changing characteristics of the object.

In the step (5), the NR-Kmeans algorithm is used to classify the elements _yi of the variable feature matrix Y, and the sparse data disturbed by noise is removed from the variable feature matrix Y to obtain a dense point set Y'; the steps are as follows:

的元素点

Q_b(y_α)表示包含元素点y_α在内的y_α的最近邻元素点集合；其中，

表示元素点y_α、

的距离，b为y_α在邻域范围内最近邻元素点个数，f＝1，2，...，b；(5.1) In the matrix Y=[y ₁ , y ₂ , ..., y _N ], for each element point y _α , y _α ∈ Y, given the neighborhood radius δ, get the

element point of

Q _b (y _α ) represents the nearest neighbor element point set of y _α including element point y _α ; where,

represents the element point y _α ,

distance, b is the number of nearest neighbor element points of y _α in the neighborhood, f=1, 2, ..., b;

(5.2) Calculate the density function value of the element point y _α; the density function value of the element point yx represents the sum of the functions of all the nearest neighbor element points in the range of the neighborhood radius δ; the calculation formula of the density function value of y _α is as follows:

为元素点y_α、

的距离，f＝1，2，...，b；In the formula,

is the element point y _α ,

distance, f = 1, 2, ..., b;

(5.3) Repeat steps (5.1) to (5.2) until the calculation of the density function values of all element points in the matrix Y is completed;

(5.4) When the density function value of y _β is greater than or equal to the average density value of all element points in the matrix Y, the element point y _β is regarded as available data and put into the dense point set Y′, that is, y _β conforms to:

In the formula, DF(y _β ) represents the density function value of y _β , and DF(y _α ) is the density function value of the data point y _α ;

(5.5) Obtain all available data in the matrix Y according to step (5.4), the available data constitute a dense point set Y', and Y' is an element point set after excluding sparse data interfered by noise.

In step (6), the dense point set Y' is divided into unchanged classes and changed classes. The gray values of the objects in the difference map corresponding to the unchanged class elements are converted to 255, and the elements in the changed classes correspond to The gray value of the object in the difference map is converted to 0, and the change detection result map is output according to the gray value of the object; the method is as follows:

(6.1) From the dense point set Y′=[y′ ₁ , y′ ₂ , ..., y′ _η ], η is the number of element points in Y′, η≤N, y′ _η means that the noise is eliminated The characteristics of the object behind the disturbed data; select the element point with the largest density value as the first initial cluster center C ₁ , and calculate the difference between the eigenvalues of other element points in Y′ and the eigenvalues of the cluster center C ₁ , select the element point with the largest difference from the _{eigenvalue of} C ₁ , that is, the point with the largest absolute value of the difference as the _second initial cluster center C _{2 ;} Class centers C ₁ and C ₂ correspond to classes D ₁ and D ₂ ;

(6.2) Calculate the difference between the element point y′ _θ in the dense point set Y′ and the eigenvalues of the two cluster centers, namely

d _θ1 = |y′ _θ -c ₁ |

d _θ2 = |y′ _θ -c ₂ |

d _θ1 represents the absolute value of the difference between the element point y' _θ and the cluster center C ₁ , d _θ2 represents the absolute value of the difference between the element point y' _θ and the cluster center C ₂ , θ=1, 2, .. ., η; divide the elements in the dense point set Y' into the class corresponding to the cluster center with the smaller absolute difference value;

(6.3) Calculate the average value c′ ₁ , c′ ₂ of the corresponding element point eigenvalues in the classes D ₁ , D ₂ , namely

In the formula, c′ ₁ , c′ ₂ are the average values of the corresponding element point eigenvalues of class D ₁ and D ₂ respectively; |D ₁ | represents the number of element points in class D ₁ , and |D ₂ | represents class D The number of element points in ₂ , y′ corresponds to the eigenvalues of the element points in the class; if c′ ₁ ≠c ₁ , update the eigenvalue c ₁ of the cluster center C ₁ of class D ₁ to c′ ₁ ; if c′ 1 ≠c 1 , ₂ ≠c ₂ , update the eigenvalue c ₂ of the cluster center C ₂ of class D ₂ to c′ ₂ ;

(6.4) Repeat steps (6.2) to (6.3) until the eigenvalues of the two cluster centers no longer change;

(6.5) Divide the two classes D ₁ , D ₂ into the unchanged class and the changed class, wherein the class with a smaller average value is the unchanged class, and the class with a larger average value is the changed class; the unchanged class The gray value of the object in the difference map corresponding to the element is converted to 255, and the gray value of the object in the difference map corresponding to the element in the changed class is converted to 0;

(6.6) Output the detection result graph according to the gray value of the object, and the area where the image changes is marked in the graph.

The invention improves the traditional kmeans clustering algorithm, and has higher detection precision and high accuracy.

The above is only the preferred embodiment of the present invention, it should be pointed out that: for those skilled in the art, without departing from the principle of the present invention, several improvements and modifications can also be made, and these improvements and modifications are also It should be regarded as the protection scope of the present invention.

Claims (8)

1. A robustness remote sensing image change detection method is characterized by comprising the following steps: the method comprises the following steps:

(1) acquiring and obtaining target area T ₁ And T ₂ The remote sensing image at a moment meets the resolution requirement and contains clear texture feature information;

(2) adopting ENVI software to realize image data preprocessing and eliminate shooting errors; using eCongintion software platform to T ₁ And T ₂ Performing multi-scale segmentation on the preprocessed images at a moment, and dividing the two images into N image blocks respectively to obtain N objects respectively; separately obtain T ₁ And T ₂ Spectral features and texture features of the image at the moment, constructing a feature vector, and obtaining a change vector set M _d ，M _d ＝[m ₁ ,m ₂ ,...m _N ]Wherein m is _N Representing a variation vector corresponding to the Nth object; will M _d The middle element is according to T ₁ Or T ₂ Arranging N objects in the time image in an arrangement mode to form a difference image;

(3) obtaining a variation vector set M according to a Relief feature selection algorithm _d Feature subspace M of _l Selecting e change characteristics with the maximum weight according to the weight of the characteristics of each object, wherein the value of e is less than the total number of the characteristics contained in the object, further respectively selecting the e change characteristics from N objects, and combining a change vector set M _d Obtaining a matrix X;

(4) performing feature dimensionality reduction on the matrix X by using a PCA algorithm to obtain a transformed change feature matrix Y, wherein Y is [ Y ═ Y [ ₁ ,y ₂ ,...,y _N ]；

(5) Using the NR-Kmeans algorithm on the elements Y of the variation feature matrix Y _i Classifying, and removing the noise from the variation characteristic matrix YObtaining dense point set Y' by using sparse data of acoustic interference; the method comprises the following steps:

(5.1) in matrix Y ═ Y ₁ ,y ₂ ,...,y _N ]For each element point y _α ，y _α E.g. Y, giving neighborhood radius delta to obtain the coincidence

Element point of (2)

Q _b (y _α ) The representation contains the element point y _α Y inside _α A set of nearest neighbor element points; wherein,

representing element point y _α 、

B is y _α The number of nearest neighbor element points in a neighborhood range, wherein f is 1, 2.

(5.2) calculating element point y _α The density function value of (a); element point y _α The density function value of (2) represents the sum of functions of all nearest neighbor element points in the neighborhood radius delta range; y is _α The calculation formula of the density function value is as follows:

in the formula,

is an element point y _α 、

1,2, …, b;

(5.3) repeating the steps (5.1) - (5.2) until the density function values of all element points in the matrix Y are calculated;

(5.4) when y _β When the density function value is more than or equal to the average density value of all the element points in the matrix Y, the element point Y is added _β Put the data considered available into the dense set of points Y', i.e. Y _β The following conditions are met:

in the formula, DF (y) _β ) Denotes y _β Value of density function, DF (y) _α ) Is a data point y _α The density function value of (a);

(5.5) obtaining all available data in the matrix Y according to the step (5.4), wherein the available data form a dense point set Y ', and Y' is an element point set from which sparse data interfered by noise is removed;

(6) and dividing the dense point set Y' into an unchanged class and a changed class, converting the gray value of the object in the difference image corresponding to the unchanged class element into 255, converting the gray value of the object in the difference image corresponding to the changed class element into 0, outputting a change detection result image according to the gray value of the object, and marking an image change area in the image.

2. The robust remote sensing image change detection method according to claim 1, characterized in that: in the step (2), the ENVI software is adopted to realize the image data preprocessing and eliminate the shooting error; using eCongintion software platform to T ₁ And T ₂ Performing multi-scale segmentation on the preprocessed images at a moment, and dividing the two images into N image blocks respectively to obtain N objects respectively; separately obtain T ₁ And T ₂ The method comprises the following steps of constructing a characteristic vector by using spectral characteristics and texture characteristics of an image at a moment, and obtaining a change vector set, wherein the steps are as follows:

(2.1) separately aligning T with ENVI software ₁ And T ₂ Performing geometric correction, image registration and atmospheric correction on the images at all times to realize image data preprocessing;

(2.2) use of eCogintion software platform to T ₁ And T ₂ Performing multi-scale segmentation on the preprocessed images at a moment, and dividing the two images into N image blocks respectively to obtain N objects respectively;

(2.3) obtaining Mean-Std characteristics of the preprocessed image as spectral characteristics of the image; extracting texture features of the preprocessed image by adopting a gray level co-occurrence matrix;

(2.4) according to the spectral features and the texture features extracted from the two preprocessed images, constructing a feature vector, and obtaining a variation vector set M _d 。

3. The robust remote sensing image change detection method according to claim 2, characterized in that: step (2.3) obtaining Mean-Std characteristics of the image as spectral characteristics of the image; the method comprises the following steps:

for each object of the image, the Mean and standard deviation Std in the spectral features are calculated according to the following formula:

wherein Mean is the Mean of the gray levels of the pixels in the object, std is the standard deviation of the gray levels of the pixels in the object, A represents the number of the pixels in the object, c _p Representing the gray scale of the p-th pixel point in the object.

4. The robust remote sensing image change detection method according to claim 2, characterized in that: extracting texture features of the image by adopting a gray level co-occurrence matrix in the step (2.3); the method comprises the following steps:

(2.3.1) any pixel point (tau) in the image ₁ ,τ ₂ ) And another pixel (tau) ₁ +△ ₁ ,τ ₂ +△ ₂ ) Form point pairs of which ₁ ,△ ₂ Is an integer; setting pixelDot (tau) ₁ ,τ ₂ ) Has a gray value of f1 (τ) ₁ +△ ₁ ,τ ₂ +△ ₂ ) If the gray value of (a) is f2, the gray value of the point pair is (f1, f2), and the maximum gray level of the image is L, the combination of f1 and f2 shares L × L groups;

(2.3.2) counting the number of occurrences of each different set of (f1, f2) values in the image, and then arranging the values into a matrix, wherein the matrix element values at the matrix positions (L, L) are the number of occurrences of the two point pairs with the gray values L;

(2.3.3) normalizing the number of occurrences of each different set of (f1, f2) values to the probability of occurrence g (f1, f2) according to the total number of occurrences of (f1, f2), the normalized square matrix being called a gray level co-occurrence matrix;

(2.3.4) selecting 6 kinds of statistic provided by the gray level co-occurrence matrix as texture features, wherein the 6 kinds of statistic are contrast, entropy, energy, equality, dissimilarity and correlation;

(2.3.5) if the matrix element value at the position (z, t) of the gray level co-occurrence matrix is g (z, t), and L is the maximum gray level of the image, the texture feature is expressed as follows:

in the formula, Con represents contrast, Ent represents entropy, Ene represents energy, Hom represents equality, Dis represents dissimilarity, and Cor represents correlation;

5. the robust remote sensing image change detection method according to claim 2, characterized in that: step (2.4) obtaining a variation vector set M according to the spectral characteristics and the texture characteristics extracted from the two preprocessed images _d (ii) a The method comprises the following steps:

(2.4.1) T obtained according to step (2.3) ₁ And T ₂ Respectively constructing characteristic vectors by the spectral characteristics and the texture characteristics of the images at the moment, and calculating the difference value of the corresponding characteristic vectors of the two images; difference value M of jth characteristic vector of ith object corresponding to two images _d (i, j), as follows:

M _d (i,j)＝M _d,1 (i,j)-M _d,2 (i,j)

in the formula, M _d,1 (i, j) is T ₁ J-th feature vector, M, of i-th object of time image _d,2 (i, j) is T ₂ The jth feature vector of the ith object of the moment image; i is less than or equal to N, wherein N is the total number of the objects in the image; j is less than or equal to S _feature ,S _feature Total number of features contained for the object; s _feature ＝S _band ×8,S _band The total wave band number of the object;

according to the difference M _d (i, j) obtaining a change vector m corresponding to the ith object _i Expressed as follows:

m _i ＝(M _d (i,1),M _d (i,2),...,M _d (i,S _feature ))

(2.4.2) repeating the step (2.4.1), calculating the difference value of the characteristic vector of each corresponding object in the two images to obtain the variation vector corresponding to each object, and forming a variation vector set M _d Expressed as follows:

M _d ＝[m ₁ ,m ₂ ,...m _N ]

will M _d The middle element is according to T ₁ And T ₂ The arrangement modes of the N objects in the time image are arranged to form a difference value image.

6. The robust remote sensing image change detection method according to any one of claims 1 to 5, characterized in that: selecting an algorithm according to the Release characteristics to obtain a variation vector set M in the step (3) _d Characteristic subspace M of _l Selecting e change characteristics with the maximum weight according to the weight of the characteristics of each object, wherein the value of e is less than the total number of the characteristics contained in the object, further respectively selecting the e change characteristics from N objects, and combining a change vector set M _d Obtaining a matrix X; the method comprises the following steps:

(3.1) vector set of variations M _d Also called feature space, extracting feature space M _d ＝[m ₁ ,m ₂ ,m ₃ ,...m _N ]Of (2) subspace M _l ＝[m ₁ ,m ₂ ,m ₃ ,...m _n ]As a training set, where N < N; the feature subspace M _l Is divided into M _l1 And M _l2 Two sets of numbers, where M _l1 Containing l1 elements, M _l2 Containing l2 elements, l1+ l2 ═ n, M _l1 ∈M _l ,M _l2 ∈M _l ；

(3.2) calculation of M _l The Euclidean distance of the change vector between every two elements in the vector form a distance matrix, namely:

where d represents the distance matrix,

represents M _l The Euclidean distance of the vector change between every two middle elements;

(3.3) in the feature subspace M _l For the change vector m of the object i _i Where i is 1,2,3, …, n, if m _i ∈M _l1 Selecting M _l1 Selecting M vectors with the minimum Euclidean distance between the vector and the vector _l2 S vectors with the largest Euclidean distance; if m is _i ∈M _l2 Selecting M _l2 Selecting M vectors with the minimum Euclidean distance between the vector and the vector _l1 S vectors with the largest Euclidean distance; namely to m _i Solving a neighbor matrix Q _i ＝[q ₁ ,q ₂ ,...,q _s ,...,q _2s ]Wherein [ q ] ₁ ,q ₂ ,...,q _s ]For s vectors with the minimum Euclidean distance, [ q ] _s+1 ,q _s+2 ,...,q _2s ]S vectors with the largest Euclidean distance, wherein s is min (l1, l 2);

(3.4) calculating the weight corresponding to each feature of the object i according to the guess neighbor correlation statistic and the guess neighbor correlation statistic, namely

w _fj ＝diffr _j -diffw _j

In the formula, w _fj Is a characteristic f _j Corresponding weight, diffr _j Is a characteristic f _j Guess the neighbor correlation statistic, diffw _j Is a characteristic f _j Guessing neighbor related statistics; diffr _j And diffw _j The calculation expression of (a) is as follows:

where, s is min (l1, l2), i denotes an object i, j denotes a jth feature, and Q is 1,2,3, …, s _i (a, j) represents [ q ] ₁ ,q ₂ ,...,q _s ]When a is s +1, s +2, …,2s, Q _i (a, j) represents [ q ] _s+1 ,q _s+2 ,...,q _2s ](ii) a When b is 1,2,3, …, s, Q _i (b, j) represents [ q ₁ ,q ₂ ,...,q _s ]When b is s +1, s +2, …,2s, Q is _i (b, j) represents [ q _s+1 ,q _s+2 ,...,q _2s ]；

(3.5) calculating the obtained feature subspace M according to the step (3.4) _l Selecting e variation characteristics with the maximum weight from the weight of the characteristic of each object, further selecting the e variation characteristics from the N objects respectively, and combining the original variation vector set M _d A matrix X is obtained, represented as follows:

in the formula, X _i ＝[x _i1 x _i2 ...x _ie ]1,2, N; matrix element x _ij And (4) representing the jth feature in the e change features of the object i, wherein the value of e is smaller than the total number of the features contained in the object.

7. The robust remote sensing image change detection method according to any one of claims 1 to 5, characterized in that: performing feature dimension reduction on the matrix X by using a PCA algorithm to obtain a transformed change feature matrix Y; the method comprises the following steps:

(4.1) calculating the mean value mu of each row of the matrix according to the matrix X obtained in the step (3) _i ,i＝1,2,...,N：

Wherein j is 1,2, E (X) _i ) Represents a mathematical expectation;

(4.2) according to the mean value μ _i The covariance matrix C of the matrix X is calculated, i.e.:

wherein, i 1,2, is, N, r 1,2, is, N, j 1,2, is, e;

(4.3) calculating an eigenvalue λ of the covariance matrix C ₁ ,λ ₂ ,...,λ _N Corresponding feature vector is v ₁ ,ν ₂ ,...,ν _N I.e. λ _i v _i ＝Cv _i ,i＝1,2,…,N；

(4.4) converting the characteristic value lambda ₁ ,λ ₂ ,...,λ _N Sorting in descending order to obtain lambda' ₁ ,λ' ₂ ,...,λ' _N The corresponding feature vector is v' ₁ ,ν' ₂ ,...,ν' _N Taking the sorted eigenvalues to form a diagonal matrix P;

(4.5) calculating the matrix X at the new feature vector v' ₁ ,ν' ₂ ,...,ν' _N And (3) obtaining a change characteristic matrix Y after dimensionality reduction by the projection, namely:

Y＝X*P＝[y ₁ ,y ₂ ,...,y _N ]

wherein, y ₁ ,y ₂ ,...,y _N Representing the changing characteristics of the object.

8. The robust remote sensing image change detection method according to any one of claims 1 to 5, characterized in that: dividing the dense point set Y' into an unchanged class and a changed class, converting the gray value of the object in the difference map corresponding to the unchanged class elements into 255, converting the gray value of the object in the difference map corresponding to the changed class elements into 0, and outputting a change detection result map according to the gray value of the object; the method comprises the following steps:

(6.1) from dense point set Y '═ Y' ₁ ,y' ₂ ,...,y' _η ]Eta is the number of element points in Y ', eta is less than or equal to N, Y' _η Representing the characteristics of the object from which the data disturbed by the noise is removed; selecting the element point with the maximum density value as a first initial clustering center C ₁ Calculating the characteristic value and the clustering center C of other element points in Y ₁ The difference of the characteristic values of (A) and (B) is selected from ₁ The element point with the largest difference of eigenvalues of (2), that is, the point with the largest absolute value of the difference, is taken as the second initial clustering center C ₂ ；C ₁ ,C ₂ Corresponding characteristic value is c ₁ ,c ₂ Cluster center C ₁ And C ₂ Corresponding class D ₁ And D ₂ ；

(6.2) calculating element points Y 'in the dense point set Y' _θ Difference from the eigenvalues of the two cluster centers, i.e.

d _θ1 ＝|y' _θ -c ₁ |

d _θ2 ＝|y' _θ -c ₂ |

d _θ1 Denotes element point y' _θ And a clustering center C ₁ Absolute value of the difference in characteristics of d _θ2 Denotes element point y' _θ And a clustering center C ₂ Is equal to 1,2, η; dividing the element points in the dense point set Y' into classes corresponding to the clustering centers with smaller absolute difference values;

(6.3) calculating class D ₁ ,D ₂ Mean value c 'of characteristic values of inner corresponding element points' ₁ ,c' ₂ I.e. by

In the formula (II), c' ₁ ,c' ₂ Are respectively of class D ₁ ,D ₂ An average value of the corresponding element point feature values; | D ₁ | represents class D ₁ Number of Medium element points, | D ₂ I represents class D ₂ The number of the middle element points, y', corresponds to the characteristic value of the element points in the class; if c' ₁ ≠c ₁ Of class D ₁ Cluster center C of ₁ Characteristic value c of ₁ Is updated to c' ₁ (ii) a If c' ₂ ≠c ₂ Class D ₂ Cluster center C of ₂ Characteristic value c of ₂ Is updated to c' ₂ ；

(6.4) repeating the steps (6.2) - (6.3) until the characteristic values of the two clustering centers are not changed any more;

(6.5) two classes D ₁ ,D ₂ The method comprises the following steps of (1) dividing the method into an unchanged class and a changed class, wherein the unchanged class is the class with a smaller average value, and the changed class is the class with a larger average value; converting the gray value of the object in the difference graph corresponding to the unchanged class elements into 255, and converting the gray value of the object in the difference graph corresponding to the changed class elements into 0;

and (6.6) outputting a detection result graph according to the gray value of the object, wherein the graph marks the image change area.

Images