CN103456018B - Remote sensing image change detection method based on fusion and PCA kernel fuzzy clustering - Google Patents

Abstract

The invention discloses a remote sensing image change detection method based on fusion and PCA kernel fuzzy clustering. The remote sensing image change detection method mainly solves the problems that in the prior art, the detection effect is not ideal, the accuracy of single-type difference image detection is low, and the application range is narrow. The method comprises the steps: (1) inputting two time phase remote sensing images X1 and X2 and conducting median filtering; (2) calculating a differential image, a logarithmic specific value image and a mean value ratio image of the two images after the filtering; (3) conducting fusion on the three images to obtain an image Xd after the fusion; (4) using a PCA method for conducting feature extraction on the images after the fusion, and obtaining a feature vector of each pixel to form a feature space matrix; (5) using a kernel-based fuzzy C mean value method for clustering the feature space matrix into two classes; (6) obtaining a final change detection result image according to the clustering result. The remote sensing change detection method has the better anti-noise performance and detection accuracy, the better effects of remote sensing images of different types can be obtained, and the remote sensing image change detection method can be applied to the field of environment monitoring and disaster evaluation.

Description

Based on the method for detecting change of remote sensing image merging with pca Kernel fuzzy clustering

Technical field

The invention belongs to technical field of image processing, it is related to method for detecting change of remote sensing image, can be used for city and plan, certainly So field such as dynamic monitoring of Disaster Assessment, Land_use change and land cover pattern.

Background technology

Remote Sensing Imagery Change Detection refers to by being compared to two width of areal different times or multiple image point Analysis, obtains the change information of atural object according to the difference between image.Remote Sensing Imagery Change Detection technology is applied successfully to Various fields, the such as mutation analysis of dynamic monitoring, forest or vegetation of environmental monitoring, Land_use change and land cover pattern, disaster are commented Estimate, agricultural investigation, cities and towns Changeement and the man-made target in military affairs monitoring and ground arms deployment analysis.

Method for detecting change of remote sensing image mainly includes 3 steps: (1) Image semantic classification；(2) construction of disparity map；(3) The extraction of change information.Wherein:

In the construction of disparity map, differential technique and ratio method as two kinds of most basic methods, respectively after to correction Two phase remote sensing images subtract each other to obtain disparity map with division operation pixel-by-pixel.Due to the taking advantage of of coherent speckle noise in sar image Property feature, ratio method than differential technique more suitable for sar image difference map construction, and ratio method is insensitive to calibration error. At present when with ratio method structural differences figure mainly in the form of logarithm ratio or average ratio.Dekker proposes using logarithm ratio Method construct disparity map,, except multiplicative noise can be converted in addition to additive noise, the excursion of also reduced value image is entered for it Go compression.Average ratio method is exactly that the neighboring mean value taking respective pixel seeks ratio again, and it has higher robustness to noise.So And this ratio method, including average ratio and logarithm ratio method exaggerate to a certain extent low gray value region change it is impossible to Truly reflect change information.If gray value is from 2 to 20 change and the change from 20 to 200, ratio method will be unable to distinguish, and The change of differential technique is just apparent.Therefore, can be existed using the construction that single ratio method or differential technique method carry out disparity map Accuracy of detection is low and the problems such as narrow application range.

The conventional method of extracting of change information is directly to difference using the clustering method such as k-means or fuzzy c-means fcm Different figure is clustered, and obtains changing class and does not change class., due to not considering space neighborhood information, its noiseproof feature is past for this method Past undesirable.

A kind of novel pca change detecting method that celik t proposes.Remote sensing for two width areal different times For image, structural differences figure first, then disparity map is divided into equal in magnitude and non-overlapping copies fritters, uses pca method Characteristic vector pickup is carried out to fritter, the neighborhood fritter then each pixel in disparity map being located is mapped in characteristic vector, Constitute new feature space, with k-means method, final change testing result figure is obtained to feature space cluster.The method letter Referred to as pca-kmeans method.This method, due to employing the data analysis technique based on piecemeal, is made with the neighborhood of each pixel For its texture, there is certain anti-noise ability.But k-means clustering method belongs to the category of hard plot, it will be every Individual sample is strictly divided into a certain class, classification clear-cut, has data or the linearly inseparable number of overlap between concrete class According to Clustering Effect undesirable.In addition pca-kmeans method is in the construction of disparity map, due to only employing single type The gray value of disparity map, the grey scale pixel value of part real change region contour and non-changing class pixel relatively, can affect The detection of region of variation, causes the problem that accuracy of detection is not high enough.

Content of the invention

Present invention aim at proposing a kind of method for detecting change of remote sensing image based on fusion and pca Kernel fuzzy clustering, Undesirable with prior art Detection results, single type disparity map accuracy of detection is low, the problem of narrow application range, preferably detects Go out region of variation.

For achieving the above object, the step of realizing of the present invention includes the following:

(1) the two width remote sensing images x not obtained in the same time are read in from areal₁And x₂, and to this two width remote sensing images Carry out 3 × 3 medium filterings；

(2) obtain medium filtering after two width images disparity map:

(2a) calculate the difference of two width image respective pixel gray values after medium filtering, and be normalized, obtain differential chart xd₁；

(2b) calculate the business of two width figure respective pixel gray values after medium filtering, and after business is taken the logarithm, then carry out normalizing Change, obtain log ratio figure xd₂；

(2c) calculate the average ratio of two width figures after medium filtering, and normalize, obtain average ratio figure xd₃；

(3) to differential chart xd₁, log ratio figure xd₂And average ratio figure xd₃Carry out 3 layers of Stationary Wavelet Transform respectively, obtain The low frequency sub-band coefficient of each width figure and high-frequency sub-band coefficient；Low frequency sub-band coefficient to this 3 width figure and high-frequency sub-band system respectively Number is merged, and carries out stationary wavelet inverse transformation to the low frequency sub-band coefficient after merging and high-frequency sub-band coefficient, is merged Image xd afterwards；

(4) with pca method, feature extraction is carried out to the image xd after merging, obtain each pixel in fused image xd Characteristic vector, constitutive characteristic space matrix q；

(5) with the Fuzzy c-means Clustering method based on core, feature space matrix q is polymerized to two different classifications a and b；

(6) correspondence according to each characteristic vector being classified as each pixel of fused image xd in feature space matrix q Relation and the cluster result of feature space matrix q, xd is divided into two corresponding with feature space matrix q cluster result not Generic a ' and b ', calculates this two different classes of averages respectively, and that larger class of average is then change class, and average is relatively That little class is then non-changing class, thus obtains final change testing result.

The present invention compared with prior art has the advantage that

1. the present invention adopts the method construct disparity map of image co-registration, and different types of remote sensing images all can be obtained preferably Testing result, solve the problems, such as that single type disparity map accuracy of detection is low, narrow application range, there is preferable robustness.

2. the feature that the present invention extracts to pca, using the fuzzy clustering method based on core, initial data is mapped to higher-dimension Feature space is clustered again, and useful feature can preferably be differentiated, extract and be amplified to the method by Nonlinear Mapping, realizes More accurately cluster, reduce further the error rate of change detection.

Brief description

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

Fig. 2 is the standards change detection figure of simulation remote sensing image data collection and this data set；

Fig. 3 is the change testing result figure with the present invention and control methods to simulation remote sensing image data collection；

Fig. 4 is the standards change detection figure of ottawa area sar image data set and this data set；

Fig. 5 is the change testing result figure with the present invention and control methods to ottawa area sar image data set；

Tu6Shi Mexico suburb remote sensing image data collection and the standards change detection figure with this data set；

Fig. 7 is the present invention and control methods change testing result figure to Mexico's suburb remote sensing image data collection.

Specific embodiment

With reference to Fig. 1, the present invention specifically comprises the following steps that

Step 1: read in the two width remote sensing images x not obtained in the same time from areal₁And x₂, and to this two width image Carry out 3 × 3 medium filterings.

1a) to x₁Each of remote sensing images pixel, chooses the rectangular window that a size is 3 × 3, by window All the gray value of pixel presses order arrangement from big to small, forms a grey-level sequence, chooses the ash in sequence centre position Angle value is as filtered output valve, thus obtaining the image x after medium filtering₁′；

1b) with same method to remote sensing images x₂Processed, obtained the image x after medium filtering₂′.

Step 2: obtain the disparity map of two width images after medium filtering.

2a) to the image x after medium filtering₁' gray value x at coordinate (i, j) place₁Figure after ' (i, j) and medium filtering As x₂' gray value x at coordinate (i, j) place₂' (i, j) carries out difference operation, obtains difference value c at coordinate (i, j) place for the matrix c (i, j):

c(i,j)=|x₁′(i,j)-x₂' (i, j) |, and then obtain difference matrix c={ c (i, j) }；

2b) normalize difference matrix, obtain differential chart xd₁Value xd at coordinate (i, j) place₁(i, j):

${\mathrm{xd}}_1(i,j)=\frac{c(i,j)-\min \left(c\right)}{\max \left(c\right)-\min \left(c\right)},$ And then obtain differential chart: xd₁={xd₁(i,j)}；

2c) to the image x after medium filtering₁' gray value x at coordinate (i, j) place₁Figure after ' (i, j) and medium filtering As x₂' gray value x at coordinate (i, j) place₂' (i, j) carries out logarithm ratio computing, obtains logarithm ratio matrix d at coordinate (i, j) place Value d (i, j):

$d(i,j)=\left|\log \frac{{x_1}^{\′}(i,j)}{{x_2}^{\′}(i,j)}\right|,$ And then obtain logarithm ratio matrix d={ d (i, j) }；

2d) normalize logarithm ratio matrix, obtain log ratio figure xd₂Value xd at coordinate (i, j) place₂(i,j)

${\mathrm{xd}}_2(i,j)=\frac{d(i,j)-\min \left(d\right)}{\max \left(d\right)-\min \left(d\right)},$ And then obtain log ratio figure: xd₂={xd₂(i,j)}；

2e) to the image x after medium filtering₁' gray value x at coordinate (i, j) place₁Figure after ' (i, j) and medium filtering As x₂' gray value x at coordinate (i, j) place₂' (i, j) carries out average ratio computing, obtains average ratio matrix e at coordinate (i, j) place Value e (i, j):

$e(i,j)=1-\min (\frac{{\mathrm{\μ}}_1(i,j)}{{\mathrm{\μ}}_2(i,j)},\frac{{\mathrm{\μ}}_2(i,j)}{{\mathrm{\μ}}_1(i,j)}),$ And then obtain average ratio matrix e={ e (i, j) },

Wherein, μ₁(i, j) and μ₂(i, j) represents image x after medium filtering respectively₁' and x₂' in centered on coordinate (i, j) 3 × 3 neighborhood windows in all pixels gray value mean value；

2f) normalization average value, than matrix, obtains average ratio figure xd₃Value xd at coordinate (i, j) place₃(i, j):

${\mathrm{xd}}_3(i,j)=\frac{e(i,j)-\min \left(e\right)}{\max \left(e\right)-\min \left(e\right)},$ And then obtain average ratio figure: xd₃={xd₃(i,j)}.

Step 3: to differential chart xd₁, log ratio figure xd₂And average ratio figure xd₃Carry out 3 layers of Stationary Wavelet Transform respectively, obtain Low frequency sub-band coefficient and high-frequency sub-band coefficient to each width figure.

3a) by differential chart xd₁Carry out 3 layers of stationary wavelet to decompose, obtain xd₁Low frequency sub-band coefficient on the 3rd decomposition layer a₁, xd₁The high-frequency sub-band coefficient in 3 directions on the 1st decomposition layerxd₁3 directions on the 2nd decomposition layer High-frequency sub-band coefficientxd₁The high-frequency sub-band coefficient in 3 directions on the 3rd decomposition layerWherein, k=lh, hl, Hh, lh represent the high-frequency sub-band of horizontal direction, and hl represents the high-frequency sub-band of vertical direction, and hh represents diagonally opposed high frequency Band；D refers to differential chart；

3b) by log ratio figure xd₂Carry out 3 layers of stationary wavelet to decompose, obtain xd₂Low frequency sub-band system on the 3rd decomposition layer Number a₂, xd₂The high-frequency sub-band coefficient in 3 directions on the 1st decomposition layerxd₂3 directions on the 2nd decomposition layer High-frequency sub-band coefficientxd₂The high-frequency sub-band coefficient in 3 directions on the 3rd decomposition layerWherein, l refers to Log ratio figure；

3c) by average ratio figure xd₃Carry out 3 layers of stationary wavelet to decompose, obtain xd₃Low frequency sub-band coefficient on the 3rd decomposition layer a₃, xd₃The high-frequency sub-band coefficient in 3 directions on the 1st decomposition layerxd₃3 directions on the 2nd decomposition layer High-frequency sub-band coefficientxd₃The high-frequency sub-band coefficient in 3 directions on the 3rd decomposition layerWherein, p refers to all Value is than figure.

Step 4: the low frequency sub-band coefficient respectively step 3 being obtained and high-frequency sub-band coefficient merge:

4a) the low frequency sub-band coefficient a to differential chart₁, the low frequency sub-band coefficient a of log ratio figure₂Low frequency with average ratio figure Sub-band coefficients a₃Merged, low frequency sub-band coefficient d ifl after being merged:

difl=a₂/4+a₃/4+a₁/2；

4b) the high-frequency sub-band coefficient of the high-frequency sub-band coefficient of log ratio figure and average ratio figure is merged, melted High-frequency sub-band coefficient d ifh after conjunction_k,n(i, j):

${\mathrm{difh}}_{k,n}(i,j)=\left\{\begin{array}{c}{\mathrm{cd}}_{k,n}^l(i,j),e_{k,n}^l(i,j)<e_{k,n}^p(i,j)\\ {\mathrm{cd}}_{k,n}^p(i,j),e_{k,n}^l(i,j)\&greaterequal;e_{k,n}^p(i,j)\end{array}\right.,$

In formula, n represents the number of plies of decomposition, n=1,2,3；Represent the high-frequency sub-band coefficient of log ratio figureValue at coordinate (i, j) place；Represent the high-frequency sub-band coefficient of average ratio figureAt coordinate (i, j) place Value；Represent the local energy coefficient of log ratio figure, it is the high-frequency sub-band coefficient of log ratio figure with coordinate The quadratic sum of all elements in 3 × 3 rectangular windows centered on (i, j)；Represent the local energy system of average ratio figure Number, it be in 3 × 3 rectangular windows centered on coordinate (i, j) for the high-frequency sub-band coefficient of average ratio figure all elements square With；

High-frequency sub-band coefficient due to differential chart contains excessive noise, so not having in the fusion of high-frequency sub-band coefficient Use it.

Step 5: stationary wavelet inverse transformation is carried out to the high-frequency sub-band coefficient after merging and low frequency sub-band coefficient, is merged Image xd afterwards；

Step 6: with pca method, feature extraction is carried out to the image xd after merging, obtain each picture in fused image xd The characteristic vector of element, constitutive characteristic space matrix q.

6a) disparity map xd after merging is divided into equal in magnitude and non-overlapping copies fritters, the size of block is h × h, h >= 2, total n fritter, each fritter is converted into column vector x^t,t=1,2,...,n；

6b) according to 4a) in n column vector obtaining, calculate mean vectorThus showing that size is h²×h²Covariance matrix: $c=\frac{1}{n}\underset{t=1}{\overset{n}{\mathrm{\&sigma;}}}\left(x^t\text{-mean}\right){(x^t-\mathrm{mean})}^t,$ Wherein, t represents transposition；

6c) Eigenvalues Decomposition is carried out to covariance matrix, obtain characteristic value and characteristic vector, by characteristic value from big to small Order arranges, and selects corresponding characteristic vector, and taking the front s column vector of characteristic vector to constitute size is h²The characteristic vector of × s Matrix v, 1≤s≤h²；

6d) h × h neighborhood fritter that each pixel in the image xd after merging is located is converted into column vector y_ε, and map To step 6c) in eigenvectors matrix v on, obtain characteristic vector v of each pixel in xd_ε=v^t(y_ε- mean), ε=1, 2 ..., h × w, wherein, mean is step 6b) in required mean vector, h is the length of fused image xd, w be merge after scheme Width as xd；

6e) according to 6d) in characteristic vector, obtain size be s × hw feature space matrix q={ v_ε, ε=1,2 ..., h×w.

Step 7: with the Fuzzy c-means Clustering method based on core, feature space matrix q is polymerized to two different classifications a And b.

7a) set clusters number as 2, initialize cluster centre c_ξ, c_ξRepresent the ξ cluster centre, ξ=1,2；

7b) set cycle counter b=0, calculate subordinated-degree matrix: u^(b)={u_ξε,

Wherein, $u_{\mathrm{\&xi;\&epsiv;}}=\frac{1/{[1-k(v_{\mathrm{\&epsiv;}},c_{\mathrm{\&xi;}})]}^{1/(m-1)}}{1/{[1-k(v_{\mathrm{\&epsiv;}},c_1)]}^{1/(m-1)}+1/{[1-k(v_{\mathrm{\&epsiv;}},c_2)]}^{1/(m-1)}}$ For the element of subordinated-degree matrix,

In formula, m is fuzzy factor, v_εIt is the ε column vector in the feature space matrix q obtaining in step 6, each row Vector is as a sample, the element u of subordinated-degree matrix_ξεIndicate the subjection degree to ξ cluster centre for the ε sample；k (v_ε,c_ξ)=exp(-||v_ε-c_ξ||²/σ²) it is gaussian kernel function, σ²> 0 it is the parameter of gaussian kernel function, exp represents fetching number；

7c) with current cluster centre and subordinated-degree matrix, each cluster centre is updated, the c after being updated ′_ξ:

$c_{\mathrm{\&xi;}}^{\&prime;}=\frac{\underset{\mathrm{\&epsiv;}=1}{\overset{h\&times;w}{\mathrm{\&sigma;}}}{u}_{\mathrm{\&xi;\&epsiv;}}^{m}k(v_{\mathrm{\&epsiv;}},c_{\mathrm{\&xi;}})v_{\mathrm{\&epsiv;}}}{\underset{\mathrm{\&epsiv;}=1}{\overset{h\&times;w}{\mathrm{\&sigma;}}}{u}_{\mathrm{\&xi;\&epsiv;}}^{m}k(v_{\mathrm{\&epsiv;}},c_{\mathrm{\&xi;}})};$

7d) with current cluster centre, subordinated-degree matrix is updated, the subordinated-degree matrix after being updated: u^(b+1)= {u_ξε}；

7e) set stopping criterion for iteration as ψ as 0.001, if max is { u^(b)-u^(b+1)< ψ then stops iteration, otherwise arranges B=b+1, c_ξ=c′_ξ, go to step 7c) and continue iteration, till meeting condition, cluster terminates, and obtains cluster centre c '₁With c′₂；

7f) by c '₁Sample for cluster centre is set to a class, by c '₂Sample for cluster centre is set to b class.

Step 8: right according to each characteristic vector being classified as each pixel of fused image xd in feature space matrix q Should be related to and feature space matrix q cluster result, xd is divided into two corresponding with feature space matrix q cluster result Different classes of a ' and b ', calculates this two different classes of averages respectively, and that larger class of average is then change class, average That class less is then non-changing class, thus obtains final change testing result.

The effect of the present invention can be further illustrated by following simulation result:

1. experiment condition

Experimental situation is: windows xp, spi, cpu pentium(r) 4, fundamental frequency 2.4ghz, software platform is matlabr2010a.

First data set is that shown in simulation remote sensing image data collection, such as Fig. 2 (a) and 2 (b), wherein Fig. 2 (a) is original Image, it is a width atm(airborne thematic mapper) the 3rd band spectrum image of image, Fig. 2 (b) is simulation Modified-image, it passes through to simulate the factor impact such as radiation characteristic of the Changes in weather of the earth and electromagnetic wave and artificially embed Region of variation obtains, and image size is 470 × 335 pixels, and gray level is 256.Fig. 2 (c) is the standards change of this data set Detection figure, including 153214 non-changing pixels and 4236 change pixels.

3rd data set be ottawa area sar image data set, and it is not being shot in the same time by two width Radarsat sar image forms, and region of variation is mainly caused by Flood Disaster Loss, shown in such as Fig. 4 (a) and 4 (b), wherein schemes 4 (a) is in May, 1997 image, and Fig. 4 (b) is the image of in August, 1997, and image size is 290 × 350 pixels, and gray level is 256.Fig. 4 (c) is the standards change detection figure of this data set, including 85451 non-changing pixels and 16049 change pixels.

3rd data set be Mexico's area remote sensing image data collection, and it is not being shot in the same time by two width Landsat7etm+ the 4th band spectrum image forms, and region of variation mainly destroys what substantial amounts of vegetation was led to by fire, As shown in Fig. 6 (a) and 6 (b), wherein Fig. 6 (a) is the image in April, 2000, and Fig. 6 (b) is the image in May, 2002, and image is big Little be 512 × 512 pixels, gray level be 256.Fig. 6 (c) is the standards change detection figure of this data set, including 236545 Non-changing pixel and 25599 change pixels.

2. experimental evaluation index

The evaluation index that experiment uses is missing inspection number, flase drop number and total error number.Wherein, missing inspection number is not detect The pixel having actually occurred change summation, flase drop number is that reality does not change but is taken as the picture that change detects The summation of element, total error number is missing inspection number and flase drop number sum.

3. experiment content and experimental result

3 image data sets are changed detect with the inventive method and existing 4 kinds of change detecting methods, this 4 kinds Control methods is respectively as follows: using difference operation structural differences figure and then is denoted as differential technique with the method for k-means cluster；Using right Number is denoted as log ratio method than computing structural differences figure and then with the method for k-means cluster；It is poor to be constructed using average ratio computing Then different figure is denoted as average ratio method with the method for k-means cluster；Celik t is in article " unsupervised change detection in satellite images using principal component analysis and k-means The method that clustering " proposes, is denoted as pca-kmeans method.

Experiment 1.

First data set is changed detect with the inventive method and existing 4 kinds of change detecting methods, parameter sets Be set to, block size h be 4, s be 3, nuclear parameter σ be 2, fuzzy factor m be 2, experimental result as shown in Figure 3 and Table 1, wherein 3 (a) It is the change testing result figure of differential technique, 3 (b) is the change testing result figure of log ratio method, 3 (c) is the change of average ratio method Change testing result figure, 3 (d) is the change testing result figure of pca-kmeans method, and 3 (e) is the change testing result figure of the present invention.

From figure 3, it can be seen that 3 (b) and 3 (c) two width in figure have more pseudo- region of variation, cause excessive flase drop Number, so that total error number is excessive；Although Fig. 3 (a) visual effect looks nice, can from the numerical statistic result of table 1 To find out, the missing inspection number of differential technique is very big, thus total error number is also larger, the method for the present invention is had less compared with differential technique Error number, and good visual effect can be obtained.

Experiment 2.

Second data set is changed detect with the inventive method and existing 4 kinds of change detecting methods, parameter sets Be set to, block size h be 3, s be 3, nuclear parameter σ be 1, fuzzy factor m be 1.4, experimental result as shown in Fig. 5 and Biao 1, wherein 5 A () is the change testing result figure of differential technique, 5 (b) is the change testing result figure of log ratio method, and 5 (c) is average ratio method Change testing result figure, 5 (d) is the change testing result figure of pca-kmeans method, the change testing result figure of 5 (e) present invention.

From fig. 5, it can be seen that missing inspection information is more in Fig. 5 (a), and noise is serious；And Fig. 5 (c) then pseudo- change information is relatively Many, flase drop is serious；The inventive method then simultaneously effective extracts change information in reduction noise.

Experiment 3.

3rd data set is changed detect with the inventive method and existing 4 kinds of change detecting methods, parameter sets Be set to, block size h be 4, s be 3, nuclear parameter σ be 1.5, fuzzy factor m be 3, experimental result as shown in Fig. 7 and Biao 1, wherein 7 A () is the change testing result figure of differential technique, 7 (b) is the change testing result figure of log ratio method, and 7 (c) is average ratio method Change testing result figure, 7 (d) is the change testing result figure of pca-kmeans method, and 7 (e) is the change testing result of the present invention Figure.

From figure 7 it can be seen that 7 (a) and 7 (c) two width in figure puppet change information are more, then missing inspection information is more for Fig. 7 (b), All can be seen that the method for the present invention can preferably suppress background information from the numerical statistic result of Fig. 7 (e) and table 1, extract Effectively change information.

With the inventive method and existing 4 kinds of change detecting methods, 3 data sets are changed with the numerical statistic detecting Result is as shown in table 1.

The numerical statistic result of table 1. change detection

As can be seen from Table 1, the flase drop number of average ratio method is many, and the missing inspection number of differential technique is many, based on single type difference The method accuracy of detection of figure is not high；The method of the present invention has less total error number than using the method for single type disparity map, Compared with pca-kmeans method, present invention incorporates the advantage of image co-registration and Kernel fuzzy clustering, decrease error number, improve Change accuracy of detection；Kernel fuzzy clustering method can preferably be differentiated, be extracted by Nonlinear Mapping and amplify useful spy Levy, different data structures effectively can be clustered, better than k-means cluster result, more accurately detect variation zone Domain.

Claims (5)

1. a kind of method for detecting change of remote sensing image based on fusion and pca Kernel fuzzy clustering, comprises the steps:

(1) the two width remote sensing images x not obtained in the same time are read in from areal₁And x₂, and this two width remote sensing images is carried out 3 × 3 medium filterings；

(2) obtain medium filtering after two width images disparity map:

(2a) calculate the difference of two width image respective pixel gray values after medium filtering, and be normalized, obtain differential chart xd₁；

(2b) calculate the business of two width figure respective pixel gray values after medium filtering, and after business is taken the logarithm, then be normalized, obtain To log ratio figure xd₂；

(2c) calculate the average ratio of two width figures after medium filtering, and normalize, obtain average ratio figure xd₃；

(3) to differential chart xd₁, log ratio figure xd₂And average ratio figure xd₃Carry out 3 layers of Stationary Wavelet Transform respectively, obtain each The low frequency sub-band coefficient of width figure and high-frequency sub-band coefficient；The low frequency sub-band coefficient of this 3 width figure is merged, gives up differential chart High-frequency sub-band coefficient, and to average ratio figure xd₃With log ratio figure xd₂High-frequency sub-band coefficient merged；To low after merging Frequency sub-band coefficients and high-frequency sub-band coefficient carry out stationary wavelet inverse transformation, the image xd after being merged；

(4) with pca method, feature extraction is carried out to the image xd after merging, obtain the feature of each pixel in fused image xd Vector, constitutive characteristic space matrix q；

(5) with the Fuzzy c-means Clustering method based on core, feature space matrix q is polymerized to two different classifications a and b；

(6) corresponding relation according to each characteristic vector being classified as each pixel of fused image xd in feature space matrix q And the cluster result of feature space matrix q, xd is divided into two inhomogeneities corresponding with feature space matrix q cluster result Other a' and b', calculates this two different classes of averages respectively, and that larger class of average is then change class, and average is less That class is then non-changing class, thus obtains final change testing result.

2. the low frequency sub-band coefficient to 3 width figures described in the method according to claim 1, wherein step (3) melts Close, carry out as follows:

Difl=a₂/4+a₃/4+a₁/2

In formula, difl is the low frequency sub-band coefficient after merging, a₁, a₂, a₃Represent differential chart xd respectively₁, log ratio figure xd₂With equal Value is than figure xd₃Carry out the low frequency sub-band coefficient after Stationary Wavelet Transform.

3. the high-frequency sub-band coefficient to 2 width figures described in the method according to claim 1, wherein step (3) melts Close, carry out as follows:

${\mathrm{difh}}_{k,n}\left(i,j\right)=\left\{\begin{array}{cc}{\mathrm{cd}}_{k,n}^l\left(i,j\right),& e_{k,n}^l\left(i,j\right)<e_{k,n}^p\left(i,j\right)\\ {\mathrm{cd}}_{k,n}^p\left(i,j\right),& e_{k,n}^l\left(i,j\right)\&greaterequal;e_{k,n}^p\left(i,j\right)\end{array}\right.$

In formula, n represents the number of plies of decomposition, n=1, and 2,3, l and p are respectively log ratio figure and average ratio figure；Represent Log ratio figure carry out Stationary Wavelet Transform after the value at coordinate (i, j) place for the high-frequency sub-band coefficient；Represent average Carry out the value at coordinate (i, j) place for the high-frequency sub-band coefficient after Stationary Wavelet Transform than figure, k=lh, hl, hh, lh represent level The high-frequency sub-band in direction, hl represents the high-frequency sub-band of vertical direction, and hh represents diagonally opposed high-frequency sub-band；Represent The local energy coefficient of log ratio figure, it be log ratio figure high-frequency sub-band coefficient centered on coordinate (i, j) 3 × 3 The quadratic sum of all elements in rectangular window；Represent the local energy coefficient of average ratio figure, it is the height of average ratio figure The quadratic sum of all elements in 3 × 3 rectangular windows centered on coordinate (i, j) for the frequency sub-band coefficients；difh_k(i, j) is to merge The value at coordinate (i, j) place for the high-frequency sub-band coefficient afterwards.

4. the image xd after merging is entered with pca method described in the method according to claim 1, wherein step (4) Row feature extraction, obtains the characteristic vector of each pixel in image xd, constitutive characteristic space matrix q, carries out as follows:

4a) disparity map xd after merging is divided into equal in magnitude and non-overlapping copies fritters, the size of block is h × h, h >=2, altogether There is n fritter, each fritter is converted into column vector x^t, t=1,2 ..., n；

4b) according to 4a) in n column vector obtaining, calculate mean vectorThus showing that size is h²×h² Covariance matrix:Wherein, t represents transposition；

4c) Eigenvalues Decomposition is carried out to covariance matrix, obtain characteristic value and characteristic vector, by characteristic value order from big to small Arrangement, selects corresponding characteristic vector, and taking the front s column vector of characteristic vector to constitute size is h²The eigenvectors matrix of × s V, 1≤s≤h²；

4d) h × h neighborhood fritter that each pixel in the image xd after merging is located is converted into column vector y_ε, and it is mapped to step On eigenvectors matrix v in 4c), obtain characteristic vector v of each pixel in xd_ε=v^t(y_ε- mean), ε=1,2 ..., h × w, wherein, mean is step 4b) in required mean vector, h is the length of fused image xd, and w is fused image xd Wide；

4e) according to 4d) in characteristic vector, obtain size be s × hw feature space matrix q={ v_ε, ε=1,2 ..., h × w.

5. described in the method according to claim 1, wherein step (5) with the Fuzzy c-means Clustering method based on core, Feature space matrix q is polymerized to two different classifications a and b, carries out as follows:

5a) set clusters number as 2, initialize cluster centre c_ξ, c_ξRepresent the ξ cluster centre, ξ=1,2；

5b) set cycle counter b=0, calculate subordinated-degree matrix: u^(b)={ u_ξε,

Wherein,For the element of subordinated-degree matrix,

In formula, m is fuzzy factor, v_εIt is the ε column vector in the feature space matrix q obtaining in step 4, each column vector As a sample；k(v_ε,c_ξ)=exp (- | | v_ε-c_ξ||²/σ²) it is gaussian kernel function, σ²> 0 be gaussian kernel function parameter, Exp represents fetching number；

5c) with current cluster centre and subordinated-degree matrix, each cluster centre is updated, the c' after being updated_ξ:

$c_{\mathrm{\&xi;}}^{\&prime;}=\frac{\underset{\mathrm{\&epsiv;}=1}{\overset{h\&times;w}{\&sigma;}}{u}_{\mathrm{\&xi;}\mathrm{\&epsiv;}}^{m}k(v_{\mathrm{\&epsiv;}},c_{\mathrm{\&xi;}})v_{\mathrm{\&epsiv;}}}{\underset{\mathrm{\&epsiv;}=1}{\overset{h\&times;w}{\&sigma;}}{u}_{\mathrm{\&xi;}\mathrm{\&epsiv;}}^{m}k(v_{\mathrm{\&epsiv;}},c_{\mathrm{\&xi;}})};$

5d) with current cluster centre, subordinated-degree matrix is updated, the subordinated-degree matrix after being updated: u^(b+1)= {u_ξε}；

5e) set stopping criterion for iteration as ψ as 0.001, if max is { u^(b)-u^(b+1)< ψ, then stop iteration, otherwise arranges b=b + 1, c_ξ=c'_ξ, go to step 5c) and continue iteration, till meeting condition, cluster terminates, and obtains cluster centre c '₁And c'₂；

5f) by c '₁Sample for cluster centre is set to a class, by c'₂Sample for cluster centre is set to b class.