CN102254323A - Method for carrying out change detection on remote sensing images based on treelet fusion and level set segmentation - Google Patents

Abstract

The invention discloses a method for carrying out change detection on remote sensing images based on treelet fusion and level set segmentation, and mainly solves the problem that much pseudo-change information exists in the existing change detection methods. The method is implemented through the following steps: inputting two time-phase remote sensing images, then respectively carrying out mean shift filtering on each image so as to obtain two time-phase filtered images; respectively carrying out two-dimensional stationary wavelet decomposition on the two time-phase filtered images three times under different level numbers; carrying out subtraction on wavelet coefficient matrixes of corresponding directional son-bands of the filtered images with the same decomposition level number; carrying out enhancement and two-dimensional wavelet inverse transformation reconstruction on wavelet coefficient difference matrixes in horizontal and vertical directions by using a sobel operator; and fusing the reconstruction images with different decomposition level numbers so as to obtain a final difference map by using a treelet algorithm, then carrying out level set segmentation on the difference map so as to obtain a change detection result. By using the method disclosed by the invention, the accuracy of the change detection result can be improved effectively, and the edge feature of a change area can be maintained better, therefore, the method can be applied to the fields of natural disaster analysis, land resource monitoring, and the like.

Description

Based on the Remote Sensing Imagery Change Detection that treelet merges and level set is cut apart

Technical field

The invention belongs to technical field of image processing, be a kind ofly to merge and level set is cut apart method for detecting change of remote sensing image specifically, can be used for the remote Sensing Image Analysis of numerous areas such as land resource monitoring, disaster analysis, urban development planning based on treelet.

Background technology

The Remote Sensing Imagery Change Detection technology is the important component part of remote sensing images research, and it compares analysis to the image of same place different times, obtains the atural object change information of people's needs according to the difference between the image.

In the change detecting method of remote sensing images, the simplest common method is directly the gradation of image value to be carried out difference to obtain difference image, utilizes threshold value to distinguish and changes class and non-variation class.Yet the remote sensing images between the different moment are because factors such as illumination, radiation cause the gradation of image value to have deviation, simply the gray scale of image is done to differ among the change detecting method result who obtains disparity map and can have more pseudo-change information, simultaneously, how to estimate accurately that threshold value is a bottleneck problem always.

In order to improve change-detection result's precision, relevant scholar has proposed many improved methods: the change detecting method that scholars such as Yakoub Bazi proposed in meeting article " A Variational Level-Set Method for Unsupervised Change Detection in Remote Sensing Images " based on the CV model in 2009, the author further improved theory part and increased experimental data on the basis of above-mentioned meeting article in 2010, had delivered journal of writings " Unsupervised Change Detection in Multispectral Remotely Sensed Imagery With Level Set Methods " again.This method does not need to estimate change threshold, but thereby the CV model is subjected to the accuracy that The noise influences the result easily.Scholars such as Turgay Celik propose a kind of change detecting method based on the dual-tree complex wavelet conversion in article " A Bayesian approach to unsupervised multiscale change detection in synthetic aperture radar image ".Because the down-sampling characteristic of dual-tree complex wavelet conversion, this method at first will be carried out interpolation to original image, the selection of interpolation method has certain influence to the result, in addition, respectively low frequency coefficient matrix and high frequency coefficient matrix are cut apart and tired easily integration class error, influence change-detection result's precision equally.

Summary of the invention

The objective of the invention is to overcome the deficiency of above-mentioned existing Remote Sensing Imagery Change Detection technology, proposed a kind of based on treelet merges and level set is cut apart method for detecting change of remote sensing image, to keep the edge of region of variation preferably, reduce pseudo-change information, improve the precision of change-detection.

For achieving the above object, detection method of the present invention comprises the steps:

(1) two width of cloth of the input size of registration is the multi-temporal remote sensing image I of m * n ¹And I ²Carry out average drifting filtering respectively, obtain image X after the filtering ¹And X ²

(2) to image X after the filtering ¹And X ²Carrying out 3 two-dimentional stationary wavelets respectively decomposes, and the highest decomposition number of plies that image decomposes for 3 times after each width of cloth filtering is respectively s (s=1, ... 3), once two-dimentional stationary wavelet decomposes and obtains four wavelet coefficient matrixes, promptly low frequency coefficient matrix with three represent respectively level, vertical, to the high frequency wavelet matrix of coefficients of angular direction;

(3) to image X after the filtering ¹And X ²Under the identical decomposition number of plies counterparty do to the wavelet coefficient matrix of subband poor, the low frequency wavelet coefficient difference matrix that obtains each decomposition layer s and three expression levels respectively, vertical, to the high frequency wavelet coefficient difference matrix of angular direction;

(4) utilize the sobel operator to strengthen to horizontal direction wavelet coefficient difference matrix in the step (3) and vertical direction wavelet coefficient difference matrix, keep low frequency wavelet coefficient difference matrix and constant angular direction wavelet coefficient difference matrix;

(5) use (3) medium and low frequency wavelet coefficient difference matrix, to level, vertical direction wavelet coefficient difference matrix after strengthening in angular direction wavelet coefficient difference matrix and (4), each decomposition layer s is carried out the contrary stationary wavelet conversion of two dimension, obtain the reconstructed image RIs of each decomposition layer s;

(6) the reconstructed image RIs to each decomposition layer s uses the treelet conversion to merge the disparity map D after obtaining merging;

(7) the disparity map D after merging is carried out level set and cut apart, obtain change-detection figure Z as a result.

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

(1) the present invention is owing to do difference and poor matrix reconstruct is come structural differences figure by the wavelet coefficient matrix that picture breakdown after two width of cloth filtering is obtained, thereby can alleviate the influence that the remote sensing images between the different moment cause the gradation of image value to exist deviation that the result is caused owing to factors such as illumination, radiation, thereby can reduce the pseudo-change information among the result, improve change-detection result's precision;

(2) the present invention is owing to adopt level, the vertical direction template of sobel operator that the wavelet coefficient difference matrix of level, vertical direction is strengthened, and makes that the edge feature of geometric configuration is more outstanding in the image;

(3) the present invention adopts three disparity map that obtain against stationary wavelet conversion reconstruct of treelet transfer pair to merge, this conversion can obtain level clustering tree and one group of multiple dimensioned orthogonal basis of a reflection data inner structure, the present invention utilizes this group orthogonal basis that the disparity map that reconstruct under three different decomposition numbers of plies obtains is merged, than directly the disparity map before merging being cut apart the result who obtains, it is higher that the disparity map after merging is cut apart the precision as a result that obtains;

(4) the present invention is owing to adopt the average drifting method that original image is carried out filtering, and the noise information in not only can the filtering image also has good edge retention performance;

(5) the present invention is owing to adopt two-dimentional stationary wavelet to decompose, and stationary wavelet has redundancy and translation invariance, therefore needn't carry out interpolation to original image, the error of having avoided interpolation algorithm to bring.

Description of drawings

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

Fig. 2 is 2 o'clock phase remote sensing images and the corresponding reference diagram that the present invention uses;

Fig. 3 is that the disparity map that obtains is merged in the present invention by treelet;

Fig. 4 is with the present invention and the control methods change-detection that experiment obtains to simulation remotely-sensed data collection figure as a result;

Fig. 5 is the change-detection that true remotely-sensed data collection experiment obtained with the present invention and control methods figure as a result.

Embodiment

With reference to Fig. 1, enforcement of the present invention is as follows:

Step 1 is to two of the input registration of phase and the size remote sensing images I that are m * n simultaneously not ¹And I ², shown in 2 (a) and Fig. 2 (b), carry out average drifting filtering respectively, obtain image X after the filtering ¹And X ²

(1a) to each pixel of input picture, according to following formula computation of mean values drift vector m _h(x) value:

$m_h\left(x\right)=\frac{{\mathrm{\Σ}}_{i=1}^{n}K\left(\frac{x-x_i}{h}\right)x_i}{{\mathrm{\Σ}}_{i=1}^{n}K\left(\frac{x-x_i}{h}\right)}$

In the formula, K (x) is a gaussian kernel function, and h is the bandwidth vector of gaussian kernel function K (x), and the present invention gets h={5, and 5}, x are the gray-scale value of current pixel point, x _iBe the gray-scale value of current pixel point bandwidth h scope interior pixel point, n is the total number that falls into current pixel point bandwidth h scope interior pixel point;

(1b) with the m that calculates _h(x) value is upgraded the gray-scale value of current pixel point;

If (1c) satisfy condition || m _h(x)-x||≤ε _s, ε _sBe given permissible error, then stop iteration, m _h(x) value is as the gray-scale value output of current pixel point, m at this moment _h(x) value is through this gray values of pixel points after the average drifting filtering, otherwise returns step (1a);

(1d) all pixels in the input picture are all carried out the operation of step (1a)-(1c), image after the output filtering.

Step 2 is to image X after the filtering ¹And X ²Carrying out 3 two-dimentional stationary wavelets respectively decomposes, and the highest decomposition number of plies that image decomposes for 3 times after each width of cloth filtering is respectively s (s=1, ... 3), once two-dimentional stationary wavelet decomposition obtains 4 wavelet coefficient matrixes, promptly 1 low frequency coefficient matrix with 3 represent respectively level, vertical, to the high frequency wavelet matrix of coefficients of angular direction, then 3 of each width of cloth image decomposition all can obtain 12 wavelet coefficient matrixes, comprise 3 low frequency coefficient matrixes and 9 high frequency coefficient matrixes.

Step 3 is to image X after the filtering ¹And X ²Under the identical decomposition number of plies counterparty do to the wavelet coefficient matrix of subband poor, the low frequency wavelet coefficient difference matrix that obtains each decomposition layer s and three expression levels respectively, vertical, to the high frequency wavelet coefficient difference matrix of angular direction, to image X after the filtering ¹And X ²Under the identical decomposition number of plies counterparty do to the wavelet coefficient matrix of subband poor, comprise to the low frequency wavelet matrix of coefficients do the difference and the high frequency wavelet matrix of coefficients is done poor, promptly

$D_s^L=|L_s^1-L_s^2|$

$D_{\mathrm{\θ},s}^H={|H}_{\mathrm{\θ},s}^1-H_{\mathrm{\θ},s}^2|$

In the formula

With

Represent image X after the filtering respectively ¹And X ²The low frequency wavelet matrix of coefficients that the s decomposition layer obtains, For

With

Low frequency wavelet coefficient difference matrix;

With

Represent image X after the filtering respectively ¹And X ²The high frequency wavelet matrix of coefficients that the s decomposition layer obtains, θ value are 0 °, 45 ° and 90 °, represent level, diagonal sum vertical direction respectively,

For

With

High frequency wavelet coefficient difference matrix.

Step 4 utilizes the sobel operator to strengthen to horizontal direction wavelet coefficient difference matrix and vertical direction wavelet coefficient difference matrix, keeps low frequency wavelet coefficient difference matrix and constant to angular direction wavelet coefficient difference matrix.

(4a) adopt the horizontal direction template of sobel operator and the high frequency wavelet coefficient difference matrix of horizontal direction

Carry out convolution, the horizontal direction wavelet coefficient difference matrix after being enhanced

(4b) adopt the vertical direction template of sobel operator and the high frequency wavelet coefficient difference matrix of vertical direction

Carry out convolution, the vertical direction wavelet coefficient difference matrix after being enhanced

Step 5 to each decomposition layer s, adopts step (3) medium and low frequency wavelet coefficient difference matrix

To angular direction wavelet coefficient difference matrix

And level, vertical direction wavelet coefficient difference matrix after strengthening in the step (4)

With

These four wavelet coefficient matrixes carry out the contrary stationary wavelet conversion of two dimension, obtain three image RIs after the reconstruct;

Step 6 merges the reconstructed image RIs use treelet conversion of each decomposition layer s, and the disparity map D after obtaining merging specifically carries out as follows;

(6a) the image RIs after the reconstruct all is transformed to the column vector RI1 ' of m * n size, RI2 ', RI3 ' forms initial sample X=[RI1 ', RI2 ', RI3 '];

(6b) number of plies of the cluster successively l=0 of definition treelet conversion, 1 ... L-1, L are the number of column vector among the initial sample X, L=3, and at the 0th layer, each variable adopts the column vector among the initial sample X to represent, initialization basis matrix B ⁰For the unit matrix of L * L and and variable indexed set Ω=1,2 ... L}, calculate the initial covariance matrix C of sample X ⁰With initial correlation matrix M ⁰, computing formula is as follows:

C _ij＝E[(s _p-Es _p)(s _q-Es _q)]

$M_{\mathrm{ij}}=\frac{C_{\mathrm{ij}}}{\sqrt{C_{\mathrm{ii}}{C}_{\mathrm{jj}}}}$

C in the formula _IjRepresent initial covariance matrix C ⁰The calculated value of the capable j row of i, s _p, s _qTwo different column vectors among the expression sample X, E represents to ask mathematical expectation, M _IjRepresent initial correlation matrix M ⁰The calculated value of the capable j row of i;

(6c) when cluster number of plies l successively ≠ 0, seek correlation matrix M ^lIn two maximum values, the correspondence position sequence number of maximal value and second largest value is designated as α and β respectively:

$(\mathrm{\α},\mathrm{\β})=\arg \underset{i,j\∈\mathrm{\Ω}}{\max }{M}_{\mathrm{ij}}^{l-1}$

Here i＜j represents correlation matrix M respectively ^lIn the row and column of arbitrary value, and only with variable indexed set Ω in carry out;

(6d) calculate Jacobi rotation matrix J:

$J=\left[\begin{array}{ccccccc}1& L& 0& L& 0& L& 0\\ M& O& M& & M& & M\\ 0& L& c_n& L& s_n& L& 0\\ M& & M& O& M& & M\\ 0& L& -s_n& L& c_n& L& 0\\ M& & M& & M& O& M\\ 0& L& 0& L& 0& L& 1\end{array}\right]$

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

Calculate J ^TThe transposition of expression rotation matrix J:

${\mathrm{\θ}}_l={\tan }^{-1}\left[\frac{C_{\mathrm{\α\α}}^l-C_{\mathrm{\β\β}}^l\±\sqrt{{(C_{\mathrm{\α\β}}^l-C_{\mathrm{\β\β}}^l)}^2+4{\left(C_{\mathrm{\α\β}}^l\right)}^2}}{2C_{\mathrm{\α\β}}^l}\right]$ And $\left|{\mathrm{\θ}}_l\right|\≤\frac{\mathrm{\π}}{4};$

In the formula,

Expression covariance matrix C ^lMiddle row, column position number is the element value of α,

Expression covariance matrix C ^lMiddle row, column position number is the element value of β,

Expression covariance matrix C ^lMiddle line position sequence number is the element value that α, column position sequence number are β;

(6e) by the orthogonal basis B after the matrix J renewal decorrelation ^l=B ^L-1J, covariance matrix C ^l=J ^TC ^L-1J, wherein B ^L-1Be the orthogonal basis before upgrading, B ^lBe the orthogonal basis after upgrading, C ^L-1Be the covariance matrix before upgrading, C ^lBe the covariance matrix after upgrading, and utilize C ^lUpgrade correlation matrix M ^L-1Be M ^l, M ^L-1Be the correlation matrix before upgrading, M ^lBe the correlation matrix after upgrading:

$M_{\mathrm{ij}}^l=\frac{C_{\mathrm{ij}}^l}{\sqrt{C_{\mathrm{ii}}^l{C}_{\mathrm{jj}}^l}}$

Wherein For upgrading back correlation matrix M ^lIn the updating value of the capable j of i row,

For upgrading back covariance matrix C ^lIn the value of the capable j of i row,

For upgrading back covariance matrix C ^lIn the value of the capable i of i row, For upgrading back covariance matrix C ^lIn the value of the capable j of j row;

(6f) back orthogonal basis matrix B is upgraded in definition ^lIn position number be that α and position number are that two column vectors of β are respectively scaling function φ _lWith Detailfunction ψ _l, the scaling vector that defines current l layer is gathered { φ _lBe scaling function φ _lScaling vector set { φ with last layer _L-1Intersection, with the position number β of difference variable from the indexed set Ω of variable remove, promptly Ω=Ω { β };

(6g) repeating step (6c) to (6f) can get orthogonal basis matrix B=[φ finally until the top l=L-1 of cluster _L-1, ψ ₁... ψ _L-1], φ wherein _L-1∈ { φ _L-1, ψ ₁Be the Detailfunction that ground floor treelet conversion obtains, ψ _L-1The Detailfunction that obtains for top treelet conversion;

(6h) initial sample X is carried out projection, i.e. P=X * B along the direction of orthogonal basis matrix B transposition ^T, and become projection vector P the image of m * n size again, the disparity map D after obtaining merging, as shown in Figure 3.

Step 7 is carried out level set to disparity map D after merging and is cut apart, and obtains change-detection Z as a result, specifically carries out as follows:

(7a) the initialization level set function is ξ ₀:

${\mathrm{\ξ}}_0(x,y)=\left\{\begin{array}{cc}-\mathrm{\ρ}& (x_a,y_a)\∈{\mathrm{\Λ}}_0-\∂{\mathrm{\Λ}}_0\\ 0& (x_a,y_a)\∈\∂{\mathrm{\Λ}}_0\\ \mathrm{\ρ}& \mathrm{\Λ}-{\mathrm{\Λ}}_0\end{array}\right.$

Λ represents the entire image territory of disparity map D to be split, Λ in the formula ₀Be any one subregion among the Λ,

Be all borderline points of this subregion, ρ is a constant, gets ρ=6 among the present invention, (x _a, y _a) coordinate of pixel arbitrarily in the presentation video;

(7b) to initial level set function ξ ₀Develop, make total energy function ε (ξ) constantly reduce when reaching the iterations 15 of regulation, to stop, this moment, the zero level collection of level set function ξ correspondence was exactly the border of region of variation, grey scale pixel value assignment 1 with region of variation, the grey scale pixel value assignment 0 of non-region of variation, the binary map that obtains is exactly final change-detection figure Z as a result, and shown in Fig. 4 (c), the evolution formula is as follows:

$\frac{\∂\mathrm{\ξ}}{\∂t}=\mathrm{\μ}[\mathrm{\Δ\ξ}-\mathrm{div}\left(\frac{\▿\mathrm{\ξ}}{|\▿\mathrm{\ξ}|}\right)]+\mathrm{\λ\δ}\left(\mathrm{\ξ}\right)\mathrm{div}\left(g\frac{\▿\mathrm{\ξ}}{|\▿\mathrm{\ξ}|}\right)+\mathrm{vg\δ}\left(\mathrm{\ξ}\right)$

λ, μ, v are constant in the formula, get λ=5 among the present invention, μ=0.2, and v=3, δ are Dirac function, Δ is represented the Laplace's operation symbol,

Expression asks gradient, div to represent to ask divergence, and g is the Boundary Detection function,

$g=\frac{1}{1+{|\▿G_{\mathrm{\σ}}*D|}^2}$

G in the formula _σThe expression standard deviation is the σ gaussian kernel, and D is for merging the back disparity map;

Total energy function ε (ξ) is expressed as internal energy function and external energy function sum, is shown below:

ε(ξ)＝μP(ξ)+ε _g，λ，v(ξ)

P in the formula (ξ) is the internal energy function, ε _{G, λ, v}(ξ) be the external energy function,

$P\left(\mathrm{\ξ}\right)={\∫}_{\mathrm{\Λ}}\frac{1}{2}{\left(\right|\▿\mathrm{\ξ}|-1)}^{2}d{x}_{a}d{y}_a$

ε _g，λ，v(ξ)＝λL _g(ξ)+vA _g(ξ)

Here L _g(ξ) be the zero level collection length of curve calculating function of level set function ξ correspondence, A _g(ξ) quicken function for curve evolvement:

$L_g\left(\mathrm{\ξ}\right)={\∫}_{\mathrm{\Λ}}\mathrm{g\δ}\left(\mathrm{\ξ}\right)|\▿\mathrm{\ξ}|d{x}_{a}d{y}_a$

A _g(ξ)＝∫ _ΛgH(-ξ)dx _ady _a

H represents step function in the formula.

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

1. experimental data

Experimental data of the present invention is one group of simulated data collection and one group of true multi-temporal remote sensing data set, wherein the original image of simulated data collection is ATM (Airborne Thematic Mapper) the 3rd band image that is positioned at Britain Feltwell village farming district of 470 * 335 pixel sizes, registration error is about 1.5 pixels, the analog variation image is to obtain by 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 remotely-sensed data collection is 2 o'clock phase Landsat5 TM+5 band images of Italian Sardinia, size is 412 * 300,256 gray levels.

2. experimental technique

The change detecting method that scholars such as method 1:Turgay Celik propose in article " A Bayesian approach to unsupervised multiscale change detection in synthetic aperture radar image " based on the dual-tree complex wavelet conversion;

The change detecting method that scholars such as method 2:Yakoub Bazi proposed in article " A Variational Level-Set Method for Unsupervised Change Detection in Remote Sensing Images " in 2010 based on the CV model;

Method 3: the inventive method.

3. experiment content and interpretation of result

Experiment 1, with distinct methods to as Fig. 2 (a) and two width of cloth shown in Fig. 2 (b) not simultaneously mutually simulated data collection carry out change-detection, the change-detection reference diagram is shown in Fig. 2 (c), the simulation remotely-sensed data collection of Fig. 2 (a) and Fig. 2 (b) is carried out result that change-detection obtains as shown in Figure 4, the change-detection that obtains for existing method 1 of Fig. 4 (a) figure as a result wherein, the change-detection that Fig. 4 (b) obtains for existing method 2 is figure as a result, the change-detection that Fig. 4 (c) obtains for the inventive method is figure as a result, as can be seen from Figure 4, there is more pseudo-change information in existing method 1 among the figure as a result, there is tangible omission in existing method 2 among the figure as a result, and the inventive method is then relatively more approaching with reference to figure 2 (c).

Experiment 2, with distinct methods to as Fig. 2 (d) and two width of cloth shown in Fig. 2 (e) not simultaneously mutually true remotely-sensed data collection carry out change-detection, the change-detection reference diagram is shown in Fig. 2 (f), the simulation remotely-sensed data collection of Fig. 2 (d) and Fig. 2 (e) is carried out result that change-detection obtains as shown in Figure 5, the change-detection that Fig. 5 (a) obtains for existing method 1 is figure as a result, the change-detection that Fig. 5 (b) obtains for existing method 2 is figure as a result, the change-detection that Fig. 5 (c) obtains for the inventive method is figure as a result, as can be seen from Figure 5, all there is many obvious noise zone, the figure as a result of the inventive method and the most approaching among the figure as a result of existing method 1 and existing method 2 with reference to figure 2 (f).

The result carries out quantitative test to the experiment change-detection, promptly compares change-detection figure and canonical reference figure as a result, calculates performance index such as total wrong number, empty scape number, omission number, and the result is as shown in table 1.

Table 1. liang group remote sensing image data change-detection evaluation of result index

From table 1, can find out more intuitively, the region of variation of the acquisition remote sensing images that the inventive method can be effectively stably, the change-detection result's that compares with other two kinds of existing methods precision is the highest.

Claims (4)

1. one kind merges based on treelet and method for detecting change of remote sensing image that level set is cut apart, comprises the steps:

(1) two width of cloth of the input size of registration is the multi-temporal remote sensing image I of m * n ¹And I ²Carry out average drifting filtering respectively, obtain image X after the filtering ¹And X ²

(2) to image X after the filtering ¹And X ²Carrying out 3 two-dimentional stationary wavelets respectively decomposes, and the highest decomposition number of plies that image decomposes for 3 times after each width of cloth filtering is respectively s (s=1, ... 3), once two-dimentional stationary wavelet decomposes and obtains four wavelet coefficient matrixes, promptly low frequency coefficient matrix with three represent respectively level, vertical, to the high frequency wavelet matrix of coefficients of angular direction;

(3) to image X after the filtering ¹And X ²Under the identical decomposition number of plies counterparty do to the wavelet coefficient matrix of subband poor, the low frequency wavelet coefficient difference matrix that obtains each decomposition layer s and three expression levels respectively, vertical, to the high frequency wavelet coefficient difference matrix of angular direction;

(4) utilize the sobel operator to strengthen to horizontal direction wavelet coefficient difference matrix in the step (3) and vertical direction wavelet coefficient difference matrix, keep low frequency wavelet coefficient difference matrix and constant angular direction wavelet coefficient difference matrix;

(5) use (3) medium and low frequency wavelet coefficient difference matrix, to level, vertical direction wavelet coefficient difference matrix after strengthening in angular direction wavelet coefficient difference matrix and (4), each decomposition layer s is carried out the contrary stationary wavelet conversion of two dimension, obtain the reconstructed image RIs of each decomposition layer s;

(6) the reconstructed image RIs to each decomposition layer s uses the treelet conversion to merge the disparity map D after obtaining merging;

(7) the disparity map D after merging is carried out level set and cut apart, obtain change-detection figure Z as a result.

2. method for detecting change of remote sensing image according to claim 1, wherein step (3) is described to image X after the filtering ¹And X ²Under the identical decomposition number of plies counterparty do to the wavelet coefficient matrix of subband poor, comprise to the low frequency wavelet matrix of coefficients do the difference and the high frequency wavelet matrix of coefficients is done poor, promptly

$D_s^L=|L_s^1-L_s^2|$

$D_{\mathrm{\θ},s}^H={|H}_{\mathrm{\θ},s}^1-H_{\mathrm{\θ},s}^2|$

In the formula

With

Represent image X after the filtering respectively ¹And X ²The low frequency wavelet matrix of coefficients that the s decomposition layer obtains,

For With

Low frequency wavelet coefficient difference matrix;

With

Represent image X after the filtering respectively ¹And X ²The high frequency wavelet matrix of coefficients that the s decomposition layer obtains, θ value are 0 °, 45 ° and 90 °, represent level, diagonal sum vertical direction respectively,

For

With

High frequency wavelet coefficient difference matrix.

3. method for detecting change of remote sensing image according to claim 1, wherein step (4) is described utilizes the sobel operator to strengthen to horizontal direction wavelet coefficient difference matrix and vertical direction wavelet coefficient difference matrix, specifically carries out as follows:

Adopt the horizontal direction template of sobel operator and the high frequency wavelet coefficient difference matrix of horizontal direction

Carry out convolution, the horizontal direction wavelet coefficient difference matrix after being enhanced

Adopt the vertical direction template of sobel operator and the high frequency wavelet coefficient difference matrix of vertical direction

Carry out convolution, the vertical direction wavelet coefficient difference matrix after being enhanced

4. method for detecting change of remote sensing image according to claim 1, wherein the image RIs of the described employing of step (6) treelet algorithm after to reconstruct merges, and carries out as follows:

(6a) the image RIs after the reconstruct all is transformed to the column vector RI1 ' of m * n size, RI2 ', RI3 ' forms initial sample X=[RI1 ', RI2 ', RI3 '];

(6b) number of plies of the cluster successively l=0 of definition treelet conversion, 1 ... L-1, L are the number of column vector among the initial sample X, L=3, and at the 0th layer, each variable adopts the column vector among the initial sample X to represent, initialization basis matrix B ⁰For the unit matrix of L * L and and variable indexed set Ω=1,2 ... L}, calculate the initial covariance matrix C of sample X ⁰With initial correlation matrix M ⁰, computing formula is as follows:

C _ij＝E[(s _p-Es _p)(s _q-Es _q)]

$M_{\mathrm{ij}}=\frac{C_{\mathrm{ij}}}{\sqrt{C_{\mathrm{ii}}{C}_{\mathrm{jj}}}}$

C in the formula _IjRepresent initial covariance matrix C ⁰The calculated value of the capable j row of i, s _p, s _qTwo different column vectors among the expression sample X, E represents to ask mathematical expectation, M _IjRepresent initial correlation matrix M ⁰The calculated value of the capable j row of i;

(6c) when cluster number of plies l successively ≠ 0, seek correlation matrix M ^lIn two maximum values, the correspondence position sequence number of maximal value and second largest value is designated as α and β respectively:

$(\mathrm{\α},\mathrm{\β})=\arg \underset{i,j\∈\mathrm{\Ω}}{\max }{M}_{\mathrm{ij}}^{l-1}$

Here i＜j represents correlation matrix M respectively ^lIn the row and column of arbitrary value, and only with variable indexed set Ω in carry out;

(6d) calculate Jacobi rotation matrix J:

$J=\left[\begin{array}{ccccccc}1& L& 0& L& 0& L& 0\\ M& O& M& & M& & M\\ 0& L& c_n& L& s_n& L& 0\\ M& & M& O& M& & M\\ 0& L& -s_n& L& c_n& L& 0\\ M& & M& & M& O& M\\ 0& L& 0& L& 0& L& 1\end{array}\right]$

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

Calculate J ^TThe transposition of expression rotation matrix J:

${\mathrm{\θ}}_l={\tan }^{-1}\left[\frac{C_{\mathrm{\α\α}}^l-C_{\mathrm{\β\β}}^l\±\sqrt{{(C_{\mathrm{\α\β}}^l-C_{\mathrm{\β\β}}^l)}^2+4{\left(C_{\mathrm{\α\β}}^l\right)}^2}}{2C_{\mathrm{\α\β}}^l}\right]$ And $\left|{\mathrm{\θ}}_l\right|\≤\frac{\mathrm{\π}}{4};$

In the formula,

Expression covariance matrix C ^lMiddle row, column position number is the element value of α,

Expression covariance matrix C ^lMiddle row, column position number is the element value of β,

Expression covariance matrix C ^lMiddle line position sequence number is the element value that α, column position sequence number are β;

(6e) by the orthogonal basis B after the matrix J renewal decorrelation ^l=B ^L-1J, covariance matrix C ^l=J ^TC ^L-1J, wherein B ^L-1Be the orthogonal basis before upgrading, B ^lBe the orthogonal basis after upgrading, C ^L-1Be the covariance matrix before upgrading, C ^lBe the covariance matrix after upgrading, and utilize C ^lUpgrade correlation matrix M ^L-1Be M ^l, M ^L-1Be the correlation matrix before upgrading, M ^lBe the correlation matrix after upgrading:

$M_{\mathrm{ij}}^l=\frac{C_{\mathrm{ij}}^l}{\sqrt{C_{\mathrm{ii}}^l{C}_{\mathrm{jj}}^l}}$

Wherein

For upgrading back correlation matrix M ^lIn the updating value of the capable j of i row,

For upgrading back covariance matrix C ^lIn the value of the capable j of i row,

For upgrading back covariance matrix C ^lIn the value of the capable i of i row,

For upgrading back covariance matrix C ^lIn the value of the capable j of j row;

(6f) back orthogonal basis matrix B is upgraded in definition ^lIn position number be that α and position number are that two column vectors of β are respectively scaling function φ _lWith Detailfunction ψ _l, the scaling vector that defines current l layer is gathered { φ _lBe scaling function φ _lScaling vector set { φ with last layer _L-1Intersection, with the position number β of difference variable from the indexed set Ω of variable remove, promptly Ω=Ω { β };

(6g) repeating step (6c) to (6f) can get orthogonal basis matrix B=[φ finally until the top l=L-1 of cluster _L-1, ψ ₁... ψ _L-1], φ wherein _L-1∈ { φ _L-1, ψ ₁The Detailfunction that obtains for ground floor treelet conversion,

ψ _L-1The Detailfunction that obtains for top treelet conversion;

(6h) initial sample X is carried out projection, i.e. P=X * B along the direction of orthogonal basis matrix B transposition ^T, and become projection vector P again the image of m * n size, the disparity map D after obtaining merging.

Images