CN102629380A - Remote sensing image change detection method based on multi-group filtering and dimension reduction - Google Patents

Abstract

The invention discloses a remote sensing image change detection method based on multi-group filtering and dimension reduction. In the prior art, precision of a change detection image and operation time are insufficient. By using the method of the invention, the above problem can be solved. The method comprises the following steps: inputting two remote sensing images with different time phases, constructing a differential difference image, acquiring a segmentation result of the differential difference image, carrying out gray scale correction on the remote sensing image so as to obtain a new differential difference image; applying a Gabor filter to obtain multi-group filtering images and forming a sample matrix; using a Treelet algorithm to carry out dimension reduction on the sample matrix, carrying out K-means (KM) on the image on which the dimension reduction is performed so as to obtain a final change detection result. Experiment shows that by using the method of the invention, change detection processing efficiency can be effectively improved, image edge information can be maintained and accuracy of the change detection can be increased. The method can be used in disaster monitoring and land utilization fields.

Description

Method for detecting change of remote sensing image based on filtering of many groups and dimensionality reduction

Technical field

The invention belongs to the digital image processing techniques field, a kind of specifically method for detecting change of remote sensing image based on filtering of many groups and dimensionality reduction is applicable to the processing and the analysis of remote sensing images.

Background technology

Change-detection based on remote sensing images is exactly from the multi-temporal remote sensing image of the same geographic area that different time obtains, the technology of analysis and definite face of land variation characteristic and process.Be widely used in fields such as agricultural, environment, city planning, national defence.Through to the not analysis of the remote sensing images of phase simultaneously of areal; Can obtain the information that this area's atural object changes; These information can be used for GIS-Geographic Information System renewal, Monitoring of Resource and Environment, target dynamic supervision and strike effect assessment etc., are the research directions of remote sensing technology.

Change detecting method to two width of cloth remote sensing images behind the registration generally is to obtain disparity map earlier, then disparity map is changed and non-variation classification.Obtain the method for disparity map because raw data is directly carried out diversity ratio, can not change data itself, information is comparatively reliable.Wherein the common classification method to region of variation and non-region of variation is the Distribution calculation classification thresholds through statistical discrepancy figure; But this method is not considered the neighborhood information of differential image, thereby is prone to violent noise section is thought that region of variation influences final accuracy of detection.The sorting technique of considering neighborhood information can overcome above shortcoming.Clustering method promptly belongs to the sorting technique of considering neighborhood information, need not set up statistical model, can effectively classify to the fuzzy part of variation and non-variation.Common clustering algorithm have fuzzy mean cluster (Fuzzy C-Means, FCM) with the K mean cluster (K-Means, KM) etc.

The FCM cluster that scholars such as Susmita Ghosh propose in article " Unsupervised Change Detection of Remotely Sensed Images using Fuzzy Clustering " is carried out the method for change-detection, can obtain classification judged result effectively to the lap of region of variation in the differential image and non-region of variation.But still the noise region that is difficult to differentiate can occur, thereby influence the result of change-detection.The PCA dimensionality reduction that utilizes that Turgay Celik scholar proposes in article " Unsupervised change detection in satellite images using principal component analysis and k-means clustering " extracts the validity feature vector of neighborhood piece in the disparity map; Thereby obtain the characteristic vector space of differential image; Utilize the K mean cluster that the characteristic vector space of differential image is classified, obtain final change-detection figure as a result according to minimum Eustachian distance at last.But the characteristic vector space of differential image needs bigger calculated amount, and image is carried out the spatial information that the piece processing can reduce region of variation, influences the result of change-detection.

Summary of the invention

The objective of the invention is to deficiency, propose a kind of method for detecting change of remote sensing image,, detect region of variation quickly and accurately to reduce calculated amount based on filtering of many groups and dimensionality reduction to above-mentioned existing method for detecting change of remote sensing image.

For realizing above-mentioned purpose, change detecting method of the present invention comprises the steps:

(1) to two width of cloth of input the corresponding grey scale pixel value in multi-temporal remote sensing image T1 and T2 locus of registration carry out difference calculating, obtain a width of cloth difference disparity map Xd;

(2) calculated difference disparity map Xd changes class and the non-classification thresholds that changes class;

(3) calculate the non-grey level histogram that changes type pixel among image T1 and the image T2 according to classification thresholds, the non-type pixel that changes among the image T1 is carried out gray correction with the non-variation class pixel among the image T2, obtain the image T3 after the image T2 gray correction; The corresponding grey scale pixel value of image T1 and image T3 locus is carried out difference calculate, obtain the difference disparity map D of a width of cloth gray correction;

(4) any with 6 frequencies any one with 5 orientation angles in the Gabor filter function parameter makes up in twos; Can obtain 30 groups of filter functions; Utilize this group filter function that the difference disparity map D of gray correction is carried out filtering and obtain 30 groups of filtered differential image groups, be designated as filtering image group matrix E ₁

(5) the difference disparity map D with gray correction joins filtering image group matrix E ₁In, obtain dimensionality reduction sample matrix E;

(6) utilize the projection matrix P of Treelet algorithm computation dimensionality reduction sample matrix E, obtain the final differential image F behind the projecting direction dimensionality reduction;

(7) utilize the K mean algorithm to carry out cluster to the final differential image F that obtains, obtain final change-detection figure Z as a result.

The present invention has the following advantages compared with prior art:

(1) the present invention carries out gray correction according to the result that presorts of difference disparity map to the non-variation type pixel of original remote sensing images, can reduce the influence of noise signal effectively, improves the accuracy of change-detection.

(2) the present invention uses many group Gabor wave filters that the difference disparity map is carried out filtering, and filtering image has dispersion and the less interior dispersion of class between bigger class, improves the cluster quality and the efficiency of algorithm of image.

(3) the present invention uses the Treelet algorithm that filtered image is carried out dimensionality reduction, obtains effective filtered, need not carry out the parameter traversal, and efficiency of algorithm will be improved.

Description of drawings

Fig. 1 is realization flow figure of the present invention;

Fig. 2 is 2 o'clock phase remotely-sensed data images and the change-detection reference diagram thereof that the present invention uses;

Fig. 3 is to the change-detection of simulation remote sensing images experiment figure as a result with the present invention;

Fig. 4 is to the change-detection of true remote sensing images experiment figure as a result with the present invention.

Embodiment

With reference to Fig. 1, implementation step of the present invention is following:

Step 1; Two width of cloth sizes of input are not mutually the remote sensing images T1 and the T2 of registration simultaneously of M * N; As shown in Figure 2, (x, the pixel gray-scale value of y) locating and

carry out difference and calculate with image T1 and image T2 space correspondence position; Obtain difference disparity map

x=1; 2 ..., M; Y=1; 2 ..., N; Wherein M is the length of remote sensing images T1, and N is the wide of remote sensing images T1.

Step 2, calculated difference disparity map Xd changes class and the non-classification thresholds T that changes class _m:

(2a) initialization classification thresholds T is the average of max pixel value and minimum pixel value among the difference disparity map Xd;

(2b) utilize initialization classification thresholds T that difference disparity map Xd is divided into two types, calculate the average m that changes type pixel value respectively _cWith the non-average m that changes type pixel value _d, with m _cAnd m _dAverage be designated as T _n, as initially type of dividing threshold value T and T _nAbsolute difference greater than convergency value T _sThe time, upgrading classification thresholds T is T _n, T _s＜＜1;

(2c) repeating step (2b) is until initialization classification thresholds T and average T _nAbsolute difference satisfy convergency value T _sThe time, obtain final classification thresholds T _m

Step 3, the non-non-variation type pixel that changes among type pixel and the image T2 among the image T1 is carried out gray correction:

(3a) according to the classification thresholds T of difference disparity map Xd _mImage T1 and T2 are divided into two types, obtain two width of cloth images variation class and non-variation class separately, calculate the non-grey level histogram p that changes type pixel among the image T1 ₁With the non-grey level histogram p that changes type pixel among the image T2 ₂

(3b) with grey level histogram p _{1 does}Reference makes grey level histogram p through the histogram coupling ₂As much as possible near grey level histogram p ₁, obtain the image T3 after the image T2 gray correction;

(3c) the corresponding grey scale pixel value of image T1 and image T3 locus is carried out difference and calculate, obtain the difference disparity map D of gray correction.

Step 4 utilizes the Gabor filter function that the difference disparity map D of gray correction is carried out filtering, obtains filtering image group matrix E ₁:

The parameter of (4a) establishing the Gabor filter function comprises six different frequency f=[0,2,4,8; 16,32] and five different directions angle φ=[0, π/3; π/6, pi/2,3 π/4]; Any with 5 orientation angles any one of 6 frequencies in the parameter is made up in twos, can obtain 30 groups of Gabor filter functions, the two-dimensional Gabor filter function can be expressed as:

$G(x,y,\mathrm{\φ},f)=\exp [-(\frac{x_{\mathrm{\φ}}^2}{{2\mathrm{\σ}}_x^2}+\frac{y_{\mathrm{\φ}}^2}{{2\mathrm{\σ}}_y^2})]\cos \left(2\mathrm{\πfx\φ}\right)$

x _φ＝xcosφ+ysinφ

y _φ＝-xsinφ+ycosφ

Wherein (x y) is the spatial domain coordinate of image, and φ is the angle of filter direction and x axle, and f is the filter center frequency, x _φ, y _φBe wave filter local coordinate, represent respectively along the direction variable of φ with perpendicular to the direction variable of φ, σ _xAnd σ _yBe respectively along x _φAnd y _φThe gaussian envelope function variance of axle utilizes these 30 groups of filter functions that the difference disparity map D of gray correction is carried out two-dimensional convolution, obtains 30 filtered images;

Be column vector (4b), obtain image sets matrix after the filtering graphical representation after each width of cloth filtering:

E ₁=[v ₁... V _Lq... V ₃₀], v wherein _LqFor the column vector of image after any width of cloth filtering is represented, lq=1 ..., 30.

Step 5 joins filtering image group matrix E with the difference disparity map D of gray correction ₁In, obtain dimensionality reduction sample matrix E:

The difference disparity map D of gray correction is represented by column vector, join image sets matrix E after the filtering ₁In, get dimensionality reduction sample matrix E:

Step 6 is utilized the projection matrix P of Treelet algorithm computation dimensionality reduction sample matrix E, obtains the final differential image F behind the projecting direction dimensionality reduction:

(6a) the covariance matrix of coefficients C of calculating dimensionality reduction sample matrix E:

$C=\left[\begin{array}{cccc}{\mathrm{<figure-callout\; id="1"\; label="C"\; filenames="HDA0000140403000000021.png"\; state="\{\{state\}\}">C</figure-callout>}}_{\mathrm{1,1}}& {\mathrm{<figure-callout\; id="1"\; label="C"\; filenames="HDA0000140403000000021.png"\; state="\{\{state\}\}">C</figure-callout>}}_{\mathrm{1,2}}& <figure-callout\; id="1"\; label="...\; C"\; filenames="HDA0000140403000000021.png"\; state="\{\{state\}\}">.</figure-callout>..& C_{}{}_{1,L}\\ C_{\mathrm{2,1}}& C_{\mathrm{2,2}}& ...& C_{2,L}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ C_{L,1}& C_{L,2}& ...& C_{L,L}\end{array}\right]---\left(1\right)$

Wherein, the Elements C among the covariance matrix of coefficients C _{I, j}Computing formula is following:

$C_{i,j}=\frac{\underset{q=1}{\overset{M\&times;N}{\mathrm{\&Sigma;}}}(v_{i,q}-\stackrel{\&OverBar;}{v_i})(v_{j,q}-\stackrel{\&OverBar;}{v_j})}{\sqrt{\underset{q=1}{\overset{M\&times;N}{\mathrm{\&Sigma;}}}{(v_{i,q}-\stackrel{\&OverBar;}{v_i})}^2\underset{k=1}{\overset{M\&times;N}{\mathrm{\&Sigma;}}}{(v_{j,q}-\stackrel{\&OverBar;}{v_j})}^2}}---\left(2\right)$

C wherein _{I, j}The calculated value of the capable j row of expression covariance matrix C i,

With Represent the average of image column vector after i and j the filtering respectively, L is the dimension of dimensionality reduction sample matrix E, L=31, and satisfied 1≤i here, j≤L;

(6b) the correlation matrix A of image sets matrix E after the calculation of filtered:

$A=\left[\begin{array}{cccc}{\mathrm{<figure-callout\; id="1"\; label="A"\; filenames="HDA0000140403000000021.png"\; state="\{\{state\}\}">A</figure-callout>}}_{\mathrm{1,1}}& {\mathrm{<figure-callout\; id="1"\; label="A"\; filenames="HDA0000140403000000021.png"\; state="\{\{state\}\}">A</figure-callout>}}_{\mathrm{1,2}}& ...& A_{1,31}\\ A_{\mathrm{2,1}}& A_{\mathrm{2,2}}& ...& A_{2,L}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ A_{L,1}& A_{L,2}& ...& A_{L,L}\end{array}\right]---\left(3\right)$

Wherein, the elements A among the correlation matrix A _{I, j}Computing formula following:

$A_{i,j}=\frac{C_{i,j}}{\sqrt{C_{i,i}{C}_{j,j}}}---\left(4\right)$

A wherein _{I, j}Be the calculated value of the capable j row of i among the correlation matrix A, C _{I, j}Be the value of the capable j row of i among the covariance matrix C, C _{I, i}Be the value of the capable i row of i among the covariance matrix C, C _{J, j}Value for the capable j row of j among the covariance matrix C;

(6c) number of plies l=0 of definition decomposition, 1 ..., L-1, when l=0, initialization and variable are image sets matrix E, i.e. S after the filtering ⁽⁰⁾=E, the difference variable is D ⁽⁰⁾Be empty set, and will with the set omega of variable can with the subscript set representations be Ω=1,2 ..., 30} will differ from variable set Ф as empty set, with quadrature Dirac base

Unit matrix as L * L dimension;

(6d) when l ≠ 0, seek two maximum among correlation matrix A values, the correspondence position sequence number of maximal value and second largest value is designated as α and β respectively:

$(\mathrm{\&alpha;},\mathrm{\&beta;})=\arg \underset{i,j\&Element;\mathrm{\&Omega;}}{\max }{A}_{i,j}^{(l-1)}---\left(5\right)$

Only consider the last triangular portions of correlation matrix A this moment, and i and j represent the row and column of arbitrary value among the correlation matrix A respectively, and only with the variable set omega in carry out;

(6e) calculate Jacobi rotation matrix J:

C wherein _n=cos (θ _l), s _n=sin (θ _l); Rotation angle θ _lBy C ^(l)=J ^TC ^(l-1)J reaches Calculate:

${\mathrm{\&theta;}}_l={\mathrm{Tan}}^{-1}\left[\frac{C_{\mathrm{\&alpha;},\mathrm{\&alpha;}}-C_{\mathrm{\&beta;},\mathrm{\&beta;}}\&PlusMinus;\sqrt{{(C_{\mathrm{\&alpha;},\mathrm{\&alpha;}}-C_{\mathrm{\&beta;},\mathrm{\&beta;}})}^2+{4C}_{\mathrm{\&alpha;},\mathrm{\&beta;}}^2}}{{2C}_{\mathrm{\&alpha;},\mathrm{\&beta;}}}\right]$ And

By the orthogonal basis B after the matrix J renewal decorrelation ^(l)=B ^(l-1)J, covariance matrix C ^(l)=J ^TC ^(l-1)J, wherein B ^(l)Be the orthogonal basis after upgrading, C ^(l)Be the covariance matrix after upgrading, and utilize covariance matrix C ^(l)Upgrade correlation matrix A according to (4) formula ^(l)

(6f) definition orthogonal basis matrix B ^(l)In α row be respectively scaling function with the column vector that β is listed as

With Detailfunction ψ _l, define the scaling vector set of current l layer

It is scaling function

Scaling vector set with last layer Intersection, with β difference variable from the variable set omega remove i.e. Ω=Ω β, join among the poor variable set Ф, i.e. Ф={ β };

(6g) repeating step (6d) to (6f) obtains projection matrix

wherein

when decomposing top l=L-1

(6h) by following formula with filtering after image sets matrix E along the direction dimensionality reduction of projection matrix P, promptly

E _n＝E×P＝[E _sE _d] (6)

E wherein _nImage sets matrix behind the expression dimensionality reduction, E _sBe image array after the filtering of main energy position projection, E _dImage sets matrix after the expression edge filtering of less important energy position projection is with E _sColumn vector becomes the identical size with original difference disparity map X, promptly obtains the final differential image F behind the dimensionality reduction.

Step 7 utilizes the K mean algorithm to carry out cluster to the final differential image F that obtains, and obtains final change-detection images Z:

(7a) collection of pixels ρ={ (m, n) } of the final differential image F of definition, 1≤m≤M, 1≤n≤N, the grey scale pixel value at pixel c place is h _c, the initialization cluster centre

G=1 ... K, wherein h _Max, h _MinBe respectively behind the dimensionality reduction maximal value and the minimum value of pixel gray scale among the final differential image F, K is classification number, K=2 here;

(7b) establish when proceeding to the H time iteration, calculate any pixel c and K cluster centre C ^HDistance B (g) (c, g)=|| h _c-C ^H(g) || and g=1 ... K, if And k ∈ [1 ..., K], then define h _c∈ ω _kType, ω _kBe the current key words sorting of pixel c, the class that is about to nearest cluster centre is as classification under this point;

(7c) recomputate K the new cluster centre in classification back:

$C^{H+1}\left(g\right)=\frac{1}{N_k}\underset{{\mathrm{\&omega;}}_k}{\mathrm{\&Sigma;}}{h}_c---\left(7\right)$

N wherein _kBe ω _kPixel number in type;

(7d) work as C ^H+1(g) ≠ C ^H(g) time, repeating step (7b) is to step (7c), if C ^H+1(g)=C ^H(g), the classification results of pixel is final change-detection figure Z as a result among final differential image F this moment.

Effect of the present invention can further specify through following experimental result and analysis:

1. experimental data

Experimental data of the present invention is one group of analog variation remote sensing images and one group of true multi-temporal remote sensing image totally two picture group pictures, each group remote sensing images testing result reference diagram that all changes.The analog variation remote sensing images are 470 * 335 pixel sizes, are positioned at the image of Airborne Thematic Mapper the 3rd wave band in Britain Feltwell village farming district that the analog variation remote sensing images are to obtain through also artificial some region of variation of embedding of factor affecting such as the Changes in weather of the simulation earth and irradiation of electromagnetic waves characteristics.True remote sensing images are that size is Landsat-5Thematic Mapper the 5th band spectrum image of 300 * 412 pixels; Two width of cloth images by the Mulargia lake, Italian Sardinia in September nineteen ninety-five and in July, 1996 is regional are formed, and variation is because the lake surface water level rises caused.

2. contrast experiment and experimental evaluation index

The method of contrast that the present invention uses is described below:

Control methods 1; Be the method that FCM cluster that scholars such as Susmita Ghosh have proposed in article " Unsupervised Change Detection of Remotely Sensed Images using Fuzzy Clustering " is carried out change-detection, can obtain classification judged result effectively the lap of region of variation in the differential image and non-region of variation.

Control methods 2; It is the orthogonal characteristic vector that utilizes the neighborhood piece among the PCA calculated difference figure that Turgay Celik scholar proposes in article " Unsupervised change detection in satellite images using principal component analysis and k-means clustering "; Thereby obtain the characteristic vector space of each pixel of differential image; Utilize the K mean algorithm that each vector of the space of feature vectors of differential image is carried out 2 classification, obtain final change-detection figure as a result according to minimum Eustachian distance at last.

What test is that the change-detection result is estimated and analyzes at last.With change-detection as a result figure and reference diagram carry out subjective vision and relatively reach objective comparison, evaluation index comprises false-alarm pixel count, omission pixel count and total erroneous pixel number, and the experiment working time of each method.

3. experiment content and analysis

Experiment 1; With distinct methods to two width of cloth not simultaneously the analog variation remote sensing images of phase carry out change-detection; Two wherein original width of cloth remote sensing images are respectively shown in Fig. 2 (a) and Fig. 2 (b); The change-detection reference diagram is shown in Fig. 2 (c), and it is as shown in Figure 3 that the analog variation remote sensing images of Fig. 2 (a) and Fig. 2 (b) are carried out the result that change-detection obtains, wherein the change-detection that obtains for control methods 1 of Fig. 3 (a) figure as a result; The change-detection that Fig. 3 (b) obtains for control methods 2 is figure as a result, and Fig. 3 (c) is the change-detection that obtains with the inventive method figure as a result.As can beappreciated from fig. 3 three kinds of method maintenances on details are all relatively good, approach Fig. 2 (c) of change-detection reference diagram.

Experiment 2; With distinct methods to two width of cloth not simultaneously the true remote sensing images of phase carry out change-detection; Two wherein original width of cloth remote sensing images are respectively shown in Fig. 2 (d) and Fig. 2 (e); Its change-detection reference diagram is shown in Fig. 2 (f), and it is as shown in Figure 4 that the true remote sensing images of Fig. 2 (d) and Fig. 2 (e) are carried out the result that change-detection obtains, wherein the change-detection that obtains for control methods 1 of Fig. 4 (a) figure as a result; The change-detection that Fig. 4 (b) obtains for control methods 2 is figure as a result, and Fig. 4 (c) is the change-detection that obtains with the inventive method figure as a result.Can find out that from Fig. 4 all there is more noise region in the testing result of two kinds of control methodss, and Fig. 2 (f) of testing result figure of the present invention and change-detection reference diagram is the most approaching.

Table 1 has been listed two groups of Remote Sensing Imagery Change Detection results' evaluation index, and table 2 is comparison working time of the present invention and control methods in the experiment of two groups of Remote Sensing Imagery Change Detection.

Table 1. Remote Sensing Imagery Change Detection evaluation of result index

Can find out more intuitively that from the evaluation index of table 1 the present invention is to the validity of Remote Sensing Imagery Change Detection.

Contrast working time of table 2. Remote Sensing Imagery Change Detection

(unit: s)	Control methods 1	Control methods 2	The inventive method
				Analog image	694.42	600.88	6.33
True picture	399.34	368.67	10.04

Can find out that from table 2 the present invention can reduce working time effectively to the change-detection of remote sensing images.

Claims (4)

1. the method for detecting change of remote sensing image based on filtering of many groups and dimensionality reduction comprises the steps:

(1) to two width of cloth of input the corresponding grey scale pixel value in multi-temporal remote sensing image T1 and T2 locus of registration carry out difference calculating, obtain a width of cloth difference disparity map Xd;

(2) calculated difference disparity map Xd changes class and the non-classification thresholds that changes class;

(3) calculate the non-grey level histogram that changes type pixel among image T1 and the image T2 according to classification thresholds, the non-type pixel that changes among the image T1 is carried out gray correction with the non-variation class pixel among the image T2, obtain the image T3 after the image T2 gray correction; The corresponding grey scale pixel value of image T1 and image T3 locus is carried out difference calculate, obtain the difference disparity map D of a width of cloth gray correction;

(4) any with 6 frequencies any one with 5 orientation angles in the Gabor filter function parameter makes up in twos; Can obtain 30 groups of filter functions; Utilize this group filter function that the difference disparity map D of gray correction is carried out filtering and obtain 30 groups of filtered differential image groups, be designated as filtering image group matrix E ₁

(5) the difference disparity map D with gray correction joins filtering image group matrix E ₁In, obtain dimensionality reduction sample matrix E;

(6) utilize the projection matrix P of Treelet algorithm computation dimensionality reduction sample matrix E, obtain the final differential image F behind the projecting direction dimensionality reduction;

(7) utilize the K mean algorithm to carry out cluster to the final differential image F that obtains, obtain final change-detection figure Z as a result.

2. method for detecting change of remote sensing image according to claim 1, wherein the said calculated difference disparity map of step (2) Xd changes class and the non-classification thresholds that changes class, calculates as follows:

(2a) initialization classification thresholds T is the average of max pixel value and minimum pixel value among the difference disparity map Xd;

(2b) utilize initialization classification thresholds T that difference disparity map Xd is divided into two types, calculate the average m that changes type pixel value respectively _cWith the non-average m that changes type pixel value _d, with m _cAnd m _dAverage be designated as T _n, as initially type of dividing threshold value T and T _nAbsolute difference greater than convergency value T _sThe time, upgrading classification thresholds T is T _n, T _s＜＜1;

(2c) repeating step (2b) is until initialization classification thresholds T and average T _nAbsolute difference satisfy convergency value T _sThe time, obtain final classification thresholds T _m

3. method for detecting change of remote sensing image according to claim 1, wherein step (3) is said carries out gray correction with the non-non-variation type pixel that changes among type pixel and the image T2 among the image T1, calculates as follows:

(3a) according to the classification thresholds T of difference disparity map Xd _mWith two width of cloth the multi-temporal remote sensing image T1 and the T2 of registration be divided into two types, obtain two width of cloth images variation class and non-variation class separately, calculate the non-grey level histogram p that changes type pixel among this image T1 ₁With the non-grey level histogram p that changes type pixel among this image T2 ₂

(3b) with grey level histogram p ₁Be reference, make grey level histogram p through the histogram coupling ₂As much as possible near grey level histogram p ₁, obtain the image T3 after the image T2 gray correction;

(3c) the corresponding grey scale pixel value in image T3 locus with multi-temporal remote sensing image T1 and after proofreading and correct carries out difference calculating, obtains the difference disparity map D of gray correction.

4. method for detecting change of remote sensing image according to claim 1, the said projection matrix P that utilizes Treelet algorithm computation dimensionality reduction sample matrix E of step (5) wherein obtains the process of the final differential image F behind the projecting direction dimensionality reduction, carries out as follows:

(6a) the covariance matrix of coefficients C of calculating dimensionality reduction sample matrix E:

$C=\left[\begin{array}{cccc}{C}_{\mathrm{1,1}}& C_{\mathrm{1,2}}& ...& C_{1,L}\\ C_{\mathrm{2,1}}& C_{\mathrm{2,2}}& ...& C_{2,L}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ C_{L,1}& C_{L,2}& ...& C_{L,L}\end{array}\right]$

Wherein, the Elements C among the covariance matrix of coefficients C _{I, j}Computing formula following

$C_{i,j}=\frac{\underset{k=1}{\overset{M\&times;N}{\mathrm{\&Sigma;}}}(v_{k,i}-\stackrel{\&OverBar;}{v_i})(v_{k,j}-\stackrel{\&OverBar;}{v_j})}{\sqrt{\underset{k=1}{\overset{M\&times;N}{\mathrm{\&Sigma;}}}{(v_{k,i}-\stackrel{\&OverBar;}{v_i})}^2\underset{k=1}{\overset{M\&times;N}{\mathrm{\&Sigma;}}}{(v_{k,j}-\stackrel{\&OverBar;}{v_j})}^2}}$

C wherein _{I, j}The calculated value of the capable j row of expression covariance matrix C i,

With

Represent the average of image column vector after i and j the filtering respectively, L is the dimension of dimensionality reduction sample matrix E, L=31, and satisfied 1≤i here, j≤L;

(6b) the correlation matrix A of image sets matrix E after the calculation of filtered:

$A=\left[\begin{array}{cccc}{A}_{\mathrm{1,1}}& A_{\mathrm{1,2}}& ...& A_{1,L}\\ A_{\mathrm{2,1}}& A_{\mathrm{2,2}}& ...& A_{2,L}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ A_{L,1}& A_{L,2}& ...& A_{L,L}\end{array}\right]$

Wherein, the elements A among the correlation matrix A _{I, j}Computing formula following:

$A_{i,j}=\frac{C_{i,j}}{\sqrt{C_{i,i}{C}_{j,j}}}$

A wherein _{I, j}Be the calculated value of the capable j row of i among the correlation matrix A, C _{I, j}Be the value of the capable j row of i among the covariance matrix C, C _{I, i}Be the value of the capable i row of i among the covariance matrix C, C _{J, j}Value for the capable j row of j among the covariance matrix C;

(6c) number of plies l=0 of definition decomposition, 1 ..., L-1, when l=0, initialization and variable are image sets matrix E, i.e. S after the filtering ⁽⁰⁾=E, the difference variable is D ⁽⁰⁾Be empty set, and will with the set omega of variable can with the subscript set representations be Ω=1,2 ..., L-1} will differ from variable set Ф as empty set, with quadrature Dirac base Unit matrix as L * L dimension;

(6d) when l ≠ 0, seek two maximum among correlation matrix A values, the correspondence position sequence number of maximal value and second largest value is designated as α and β respectively:

$(\mathrm{\&alpha;},\mathrm{\&beta;})=\arg \underset{i,j\&Element;\mathrm{\&Omega;}}{\max }{A}_{i,j}^{(l-1)}$

Only consider the last triangular portions of correlation matrix A this moment, and i and j represent the row and column of arbitrary value among the correlation matrix A respectively, and only with the variable set omega in carry out;

(6e) calculate Jacobi rotation matrix J:

C wherein _n=cos (θ _l), s _n=sin (θ _l); Rotation angle θ _lBy C ^(l)=J ^TC ^(l-1)J reaches

Calculate:

${\mathrm{\&theta;}}_l={\mathrm{Tan}}^{-1}\left[\frac{C_{\mathrm{\&alpha;},\mathrm{\&alpha;}}-C_{\mathrm{\&beta;},\mathrm{\&beta;}}\&PlusMinus;\sqrt{{(C_{\mathrm{\&alpha;},\mathrm{\&alpha;}}-C_{\mathrm{\&beta;},\mathrm{\&beta;}})}^2+{4C}_{\mathrm{\&alpha;},\mathrm{\&beta;}}^2}}{{2C}_{\mathrm{\&alpha;},\mathrm{\&beta;}}}\right]$ And

By the orthogonal basis B after the matrix J renewal decorrelation ^(l)=B ^(l-1)J, covariance matrix C ^(l)=J ^TC ^(l-1)J, wherein B ^(l)Be the orthogonal basis after upgrading, C ^(l)Be the covariance matrix after upgrading, and utilize covariance matrix C ^(l)Upgrade correlation matrix A according to following formula ^(l):

$A_{i,j}^{\left(l\right)}=\frac{C_{i,j}^{\left(l\right)}}{\sqrt{C_{i,i}^{\left(l\right)}{C}_{j,j}^{\left(l\right)}}};$

(6f) definition orthogonal basis matrix B ^(l)In α row be respectively scaling function with the column vector that β is listed as

With Detailfunction ψ _l, define the scaling vector set of current l layer

It is scaling function

Scaling vector set with last layer

Intersection, with β difference variable from the variable set omega remove i.e. Ω=Ω β, join among the poor variable set Ф, i.e. Ф={ β };

(6g) repeating step (6d) to (6f) obtains projection matrix

wherein

when decomposing top l=L-1

(6h) by following formula with filtering after image sets matrix E along the direction dimensionality reduction of projection matrix P, promptly

E _n＝E×P＝[E _sE _d]

E wherein _nImage sets matrix behind the expression dimensionality reduction, E _sBe image array after the filtering of main energy position projection, E _dImage sets matrix after the expression edge filtering of less important energy position projection is with E _sColumn vector becomes the identical size with original difference disparity map X, promptly obtains the final differential image F behind the dimensionality reduction.

Images