CN103426158A - Method for detecting two-time-phase remote sensing image change - Google Patents

Abstract

The invention discloses a method for detecting two-time-phase remote sensing image change. In the method for detecting two-time-phase remote sensing image change, a gaussian process classifier is combined with a Markov random field model, effective utilization of image space information is reinforced, and change detecting precision is improved.

Description

The method of 2 o'clock phase Remote Sensing Imagery Change Detection

Technical field

The present invention relates to the Radar Technology field of remote sensing image processing, relate in particular to the method for a kind of 2 o'clock phase Remote Sensing Imagery Change Detection.

Background technology

It is one of most important application direction in field of remote sensing image processing that remote sensing images are changed to detection, and it has important using value at numerous Military and civil fields such as Natural calamity monitoring, land and resources program management, military target strike assessment.Along with the resolution of remote sensing image data increases, change and detect the change information that can obtain and become increasingly abundant, make to change the practical ranges detected and further enlarged.In prior art, remote sensing images are changed to the method detected and mainly contain: Gaussian process method and Markov random field model method.Below respectively two kinds of methods are described:

Gaussian process (Gaussian Process, be called for short GP) as a kind of core learning machine with probability meaning, for the study of kernel function provides a kind of theoretical foundation that both had, can be used for again the probability model of practice, and a set of complete theoretical frame is being provided aspect selection, study and the prediction of model.It can not only mean with the form of prior probability the priori of process, improves the process model performance, and can provide the precision parameter of the prediction of output.With support vector machine classifier, compare, the Gaussian process sorter has Bayesian formula completely and means, model parameter obviously reduces, and parameter optimization is relatively easy, more easily convergence.But the Gaussian process sorter it has been generally acknowledged that each unique point is independently, do not exist and influence each other between the further feature around unique point and its in neighborhood, this hypothesis often causes in testing result having a large amount of noises, causes accuracy of detection not high.

Markov random field model can effectively utilize the spatial context information of image.Change detecting method based on Markov random field (Markov Random Field is called for short MRF) is explained spatial context information by the form of energy function, to fall low noise impact, improves the accuracy rate of testing result.But Markov random field is more responsive to initial value, easily cause final detection output area locally optimal solution, cause the whole detection precision not high.

Summary of the invention

(1) technical matters that will solve

For solving above-mentioned one or more problems, the invention provides the method for a kind of 2 o'clock phase Remote Sensing Imagery Change Detection, to overcome unique point in the Gaussian process method and not exist and influence each other between other unique points in neighborhood around it, and the Markov random field model method to initial value than more sensitive defect, improve 2 o'clock phase Remote Sensing Imagery Change Detection accuracy of detection.

(2) technical scheme

According to an aspect of the present invention, a kind of 2 o'clock phase method for detecting change of remote sensing image are provided.The method comprises: respectively the remote sensing images of 2 o'clock phases carried out to over-segmentation, obtain a plurality of overdivided regions after dividing processing, the remote sensing images that 2 o'clock phase remote sensing images are two width the same area different times; Respectively to per a period of time the phase remote sensing images take pending overdivided region and extract color characteristic, textural characteristics as unit, and entropy feature; Color, textural characteristics, the entropy feature in corresponding pending zone in 2 o'clock phase remote sensing images are asked to poor, and required difference is as the proper vector of corresponding pending overdivided region in the error image of 2 o'clock phase remote sensing images; Utilize housebroken Gaussian process sorter, according to the characteristic of correspondence vector, overdivided region in error image is changed class and do not changed two classes classification of class, obtain the class probability of each overdivided region, and by this classification probability assignment to all pixels in this pending overdivided region on error image; Build Markov random field model on error image, initial value using the class probability value of each pixel in each overdivided region as Markov random field model, calculate the Markov random field model energy function, by energy function optimization, obtain the result of 2 o'clock phase Remote Sensing Imagery Change Detection.

(3) beneficial effect

From technique scheme, can find out, the present invention has following beneficial effect:

(1) Markov random field model is combined with the Gaussian process sorter, overcome the Gaussian process sorter and lacked the shortcoming to effective utilization of image space information, effectively utilize spatial context information and improved the variation accuracy of detection;

(2) Gaussian process and Markov random field model are combined, effectively solved the distribution form of observing in the Remote Sensing Imagery Change Detection and be difficult to the definite problem of parametrization, improved accuracy of detection.

The accompanying drawing explanation

The schematic flow sheet of phase method for detecting change of remote sensing image when Fig. 1 is the embodiment of the present invention two.

The criterion schematic diagram of in the phase method for detecting change of remote sensing image, 2 o'clock each overdivided region corresponding to phase remote sensing images being adjusted when Fig. 2 is the embodiment of the present invention two;

Phase remote sensing images when Fig. 3 is two width two, wherein, Fig. 3 a is photographic images in 2002; Fig. 3 b is photographic images in 2003;

Fig. 4 be the present invention while adopting distinct methods to two width two as shown in Figure 3 the phase remote sensing images change the result detected.Wherein, the variation testing result that Fig. 4 a is the first comparative approach, Fig. 4 b is the variation testing result of the second comparative approach, Fig. 4 c is the variation testing result of the third comparative approach, Fig. 4 d is the variation testing result that the present invention is based on the Gaussian process method for detecting change of remote sensing image that spatial context is relevant, and Fig. 4 e is artificial testing result.

Embodiment

For making the purpose, technical solutions and advantages of the present invention clearer, below in conjunction with specific embodiment, and, with reference to accompanying drawing, the present invention is described in more detail.

It should be noted that, in accompanying drawing or instructions description, similar or identical part is all used identical figure number.And in the accompanying drawings, embodiment is to simplify or convenient the sign.Moreover the implementation that does not illustrate in accompanying drawing or describe, be form known to a person of ordinary skill in the art in affiliated technical field.In addition, although this paper can provide the demonstration of the parameter that comprises particular value, should be appreciated that, parameter is without definitely equaling corresponding value, but can in acceptable error margin or design constraint, be similar to corresponding value.

In the phase method for detecting change of remote sensing image, the Gaussian process sorter is combined with Markov random field model during the present invention two, strengthen the effective utilization to image space information, improve and change accuracy of detection.

In one exemplary embodiment of the present invention, a kind of 2 o'clock phase method for detecting change of remote sensing image are provided.The schematic flow sheet of phase method for detecting change of remote sensing image when Fig. 1 is the embodiment of the present invention two.As shown in Figure 1, the present embodiment comprises the following steps:

Step S102, carry out over-segmentation to the remote sensing images of 2 o'clock phases respectively, obtains a plurality of overdivided regions after dividing processing;

Wherein, the remote sensing images of 2 o'clock above-mentioned phases are preferably the remote sensing image that original input picture obtains for the same sensor of two width the same area different times after registration process, its image spatial resolution is higher than 2 meters/pixel, and image comprises tri-wave bands of RGB.

In this step, utilize watershed segmentation method (J.Roerdink and A.Meijster, " The Watershed Transform:Definitions, Algorithms, and Parallellization Strategies; " Fundamenta Informaticae, vol.41, pp.187-228,2000.), while take each, the CIELab color characteristic of phase images is basis, and during respectively to each, the remote sensing images of phase carry out the over-segmentation processing.Each overdivided region possesses certain regional homogeneity consistance, is homogeneous region.In the over-segmentation process, in order to guarantee regional most homogeneous and to change the actual decipher detected, the chosen area area is as constraint condition, and when total number of pixels of certain cut zone is less than 200, this zone stops cutting apart;

Step S104, adjusted 2 o'clock each overdivided region corresponding to phase remote sensing images, makes the overdivided region of correspondence of 2 o'clock phase remote sensing images consistent;

In this step, for an overdivided region A in the first phase image, its corresponding region in second o'clock phase images is overdivided region B, and the adjustment criterion of overdivided region A and overdivided region B is:

If A=B, overdivided region A and overdivided region B all remain unchanged, as shown in Figure 2;

If A ≠ B, and A ∩ B=A, overdivided region A remains unchanged, and overdivided region B is divided into to B=B1 ∪ B2, wherein, the corresponding overdivided region that B1 is overdivided region A on second o'clock phase images, B2 is B-B1, as shown in Figure 2;

If A ∪ is B>A, and A ∪ B>B, overdivided region A being divided into to A-(A ∩ B) and A ∩ B two parts, overdivided region B is divided into B-(A ∩ B) and A ∩ B two parts, as shown in Figure 2.

Step S106, respectively to per a period of time the phase remote sensing images take pending overdivided region and extract color characteristic, textural characteristics as unit, and entropy feature;

This step can be divided into again step by step following:

S106a step by step, the pixel in the overdivided region target of take is unit, extracts CIELab color characteristic, Gabor textural characteristics and the entropy feature of each pixel;

Wherein, CIELab color characteristic extraction step is: original input picture is carried out to smothing filtering, each pixel p (i in image, j), centered by coordinate (i, j), get the image-region of 3*3 size, in zoning the spectrum average of pixel as the filtering of current pixel p (i, j) after spectral value, 2 o'clock phase remote sensing images after obtaining filtering and processing; By image RGB color space projection CIELab color space, obtain the color feature value of each pixel in the CIELab color space in image; The over-segmentation image that the back of take obtains is basis, and the average of the Lab color characteristic of all pixels CIELab feature regional as this in each overdivided region obtains the 3 Vc IELab color characteristics in each zone; Wherein, the RGB color space to the conversion computing formula of CIELab color space is:

$X=\frac{0.412453*R+0.357580*G+0.180423*B}{\text{0.950456}}$

$Y=\frac{0.212671*R+0.715160*G+0.072169*B}{1}$

$Z=\frac{0.019334*R+0.119193*G+0.950227*B}{1.088754}$

XT＝X＞th，YT＝Y＞th，ZT＝Z＞th，th＝0.008856

fX＝XT.*X ^1/3+(～XT).*(7.787.*X+16/116)

fY＝YT.*Y ^1/3+(～YT).*(7.787.*Y+16/116)

fZ＝ZT.*Z ^1/3+(～ZT).*(7.787.*Z+16/116)

L＝YT.*(116*Y ^1/3-16.0)+(～YT).*(903.3.*Y)

a＝500*(fX-fY)

b＝200*(fY-fZ) (1)

Wherein, X, Y, Z, fX, fY, fZ its only for calculating intermediate value, there is no concrete physical meaning.R, G and B are respectively the R/G/B component value of the RGB colour system of respective pixel.L, a and b are respectively L, a of CIELab colour system of the respective pixel after conversion and the value of each component of b.Th is judgment threshold.XT, YT, ZT mean that respectively X, Y, Z are greater than the part of judgment threshold th.

Gabor texture feature extraction step is: original input remote sensing images are converted into to gray level image, and method for transformation is for directly to average the spectral value addition of each pixel on image R, G, tri-passages of B, as the spectral value of gray level image; Structure Gabor wave filter is as follows:

x′＝xcosθ+ysinθ

y′＝-xsinθ+ycosθ (2)

Wherein, λ is wavelength, the strip direction that θ is the Gabor wave filter,

For the phase parameter of Gabor wave filter, γ is length breadth ratio, controls the index ellipsoid of Gabor wave filter.The Gabor energy function is defined as follows form:

$e_{\mathrm{\λ},\mathrm{\θ}}(x,y)=\sqrt{{\mathrm{\γ}}_{\mathrm{\λ},\mathrm{\θ},0}^2(x,y)+{\mathrm{\γ}}_{\mathrm{\λ},\mathrm{\θ},-\frac{\mathrm{\π}}{2}}^2(x,y)}---\left(3\right)$

Wherein,

With

Be respectively g _{λ, θ, 0}(x, y) and Result with the gray level image convolution; The scale parameter λ that sets the Gabor wave filter in the present invention gets respectively [12345], direction parameter θ gets respectively [0 π/4 (2* π)/4 (3* π)/4 (4* π)/4 (5* π)/4 (6* π)/4 (7* π)/4], can construct Gabor wave filter in 40, use each wave filter to gray level image filtering, be about to each Gabor wave filter and gray level image and carry out convolution; Take the image overdivided region as unit, and the average of the Gabor textural characteristics of all pixels Gabor feature regional as this in zone, obtain 40 of each zone and tie up the Gabor textural characteristics;

The entropy characteristic extraction step is: original input remote sensing images are converted into to gray level image, and method for transformation is for directly to average the spectral value addition of each pixel on image R, G, tri-passages of B, as the spectral value of gray level image; Take each zone in the over-segmentation image is unit, and the entropy in computed image zone obtains the 1 dimension entropy feature in each zone.The computing formula of entropy feature is as follows:

$\mathrm{entry}=\underset{i=1}{\overset{255}{\mathrm{\Σ}}}{p}_i\log p_i.---\left(4\right)$

Wherein, p _iThe number of pixels that in the presentation video zone, gray-scale value is i accounts for the ratio of number of pixels in whole image-region.

S106a, calculate the feature of interior each characteristic mean of pixel of overdivided region as this pending overdivided region step by step.

Step S108, ask poor by color, textural characteristics, the entropy feature in corresponding pending zone in 2 o'clock phase remote sensing images, and required difference is as the proper vector of corresponding pending overdivided region in the error image of 2 o'clock phase remote sensing images;

This step can be divided into again step by step following:

S108a, carry out normalized to CIELab color characteristic, Gabor textural characteristics, entropy feature respectively step by step, obtains normalization CIELab color characteristic, Gabor textural characteristics and entropy feature afterwards.Wherein, the normalized formula is expressed as follows:

f′＝(f-min(min(f)))/(max(max(f))-min(min(f))). (5)

Wherein, f is above-mentioned CIELab color characteristic, Gabor textural characteristics, entropy feature.

S108b step by step, according to the over-segmentation image after adjusting, on overdivided region, feature is subtracted each other one by one, and a plurality of features of obtaining pending overdivided region in 2 o'clock phase remote sensing images are poor, the proper vector of corresponding pending overdivided region in the poor formation error image of the plurality of feature.

Step S110, utilize housebroken Gaussian process sorter, according to proper vector, overdivided region in error image is changed class and do not changed two classes classification of class, obtain the class probability of each overdivided region, and by this classification probability assignment to all pixels in this pending overdivided region on error image;

This step can be divided into following steps again:

Step S110a, obtain training sample set

Wherein, X _iFor some training samples, mean the remote sensing image (image spatial resolution is higher than 2 meters/pixel, and image comprises tri-wave bands of RGB) that two width obtain through the same sensor of areal different time of registration process;

For training sample X _iCorresponding two-value classification, indication training sample X _iIn each zone belong to and change class or belong to the probability that does not change class; According to S102-S106 to each training sample X _iExtract feature, obtain about training sample X _iTake the characteristic set that zone is unit, combine each training sample data class mark Y _iComposing training sample characteristics collection

Test data set

Wherein, X _i ^*For some test sample books, the remote sensing image (image spatial resolution is higher than 2 meters/pixel, and image comprises tri-wave bands of RGB) that each test sample book is obtained through the same sensor of the same area different time of registration process by two width forms.According to S102-S106, each test sample book is extracted to feature, obtain the feature set that zone is unit of take about test sample book

Step S110b, utilize the proper vector x of the concentrated training sample of training sample _iWith class mark y _iThe Gaussian process sorter is trained, obtained the gaussian kernel parameter;

Suppose that Gauss's priori average is 0, by training sample, calculate Gauss's prior variance matrix.The present invention uses gaussian radial basis function kernel function (Gaussian Radial Basis Function, RBF), and expression formula is as follows:

$k(x_i,x_j)={\mathrm{\θ}}_0\exp (-\frac{{|x_i-x_j|}^2}{2{\mathrm{\σ}}^2})---\left(6\right)$

Wherein, θ ₀Control parameter for the amplitude of gaussian radial basis function kernel function, the yardstick that σ is the gaussian radial basis function kernel function is controlled parameter; x _i, x _iMean the characteristics of image obtained in S108a, θ ₀Be with σ two gaussian kernel parameters that need.

Step S110c, utilize Gaussian process sorter corresponding to gaussian kernel parameter, to the test feature collection that pending overdivided region proper vector forms in described error image

Classified, obtained the class probability that test feature is concentrated each pending overdivided region;

For test sample book x _i ^*(i=l+1 ... n), the present invention's acceptance of the bid remembers that an output probability value corresponding to each test sample book is y _i ^*(i=l+1 ... l+u), this probable value means the class probability of the corresponding overdivided region of this test sample book.Concrete steps are as follows:

1,, by Gauss's priori average and variance matrix, obtain implicit function f (x _i) expression formula.Suppose implicit function f (x _i) the obedience Gaussian process: f (x _i, θ)～GP (0, K), make f:x _i→ y _i(i=1,2 ..., l), wherein, y _iMean each regional corresponding class mark in training set.Therefore, can utilize the training sample feature set

Mapping relationship f by supposition: x _i→ y _i(i=1,2 ..., l), training obtains implicit function f (x _i) expression-form.

2, utilize the implicit function f (x obtained _i) expression formula, calculate the test sample book feature set

Correspondence prediction output probability y _i ^*.Wherein, prediction output probability calculation expression:

$p\left(y_i^{*}\right|x_i^{*},x_i,y_i)=\∫p\left(y_i^{*}\right|f^{*})p\left(f^{*}\right|x_i^{*},X_l,Y_l){\mathrm{df}}^{*}---\left(2\right)$

In formula

Mean noise model,

For process model, mean sample

Hidden variable f ^*Probability distribution.

3, to noise model

Due to sample in Gauss's sorter

Belong to the class mark

Probability be and implicit function f ^*Be directly proportional, so noise model

Can it be defined as:

$p(y_i^{*}=+1|f_i^{*})=\mathrm{\Φ}\left(f\left(x_i^{*}\right)\right)---\left(3\right)$

In formula, the Φ function is S type (sigmoid) function class.

4, for process model

By marginalisation f=[f ₁, f ₂..., f _l], we can obtain:

$p\left(f_i^{*}\right|x_i^{*},X_l,Y_l)=\∫p\left(f_i^{*}\right|x_i^{*},x,f)p\left(f\right|x,y)\mathrm{df}---\left(4\right)$

In formula, p (f|x, y) is the posterior probability of hidden variable.In conjunction with Gaussian process priori p (f) and likelihood function p (y|f) thereof, we can obtain posteriority and distribute:

p(f|x，y)＝p(y|f)p(f|X)/p(y|X) (5)

Due to when the given latent function f, observation data is separate, and likelihood function p (y|f) can be described as:

$p\left(y\right|f)=\underset{i=1}{\overset{n}{\mathrm{\Π}}}p\left(y_i\right|f_i)=\underset{i=1}{\overset{n}{\mathrm{\Π}}}\mathrm{\Φ}\left(y_i{f}_i\right)---\left(6\right)$

Marginal probability p (y|X) can be described as:

$p\left(y\right|X)=\∫p\left(f\right|X)\underset{i=1}{\overset{n}{\mathrm{\Π}}}p\left(y_i\right|f_i)\mathrm{df}---\left(7\right)$

5, approximate solution: due to noise model p (y in formula (6) _i| f _i) what adopt is S type (sigmoid) class function, although prior distribution p is (f _i) be Gaussian function, but posteriority distributes

And prediction distribution

Be non-Gaussian distribution form, cause Gauss's expression formula (formula (2)) of classifying to be solved by the mode of resolving.In the present invention, adopt Laplce's method of approximation to noise model

Carry out approximate solution, and the Approximate prediction distribution that distributes and provide test sample book by this approximate posteriority.

Step S110d, give all pixels in this overdivided region by the class probability assignment of each overdivided region, as changing the initial result detected.

Step S112, build Markov random field model on error image, initial value using each pixel class probability value of each overdivided region as Markov random field model, calculate the Markov random field model energy function, by energy function optimization, obtain the result of 2 o'clock phase Remote Sensing Imagery Change Detection.

This step can be divided into again step by step following:

S112a step by step, a Markov random field model G (V of structure on error image, E), in error image, each pixel is expressed as a node V in Markov random field model, the pixel of its peripheral region of each node V forms the neighborhood set N of this node V, the poor E of the limit as Markov random field model of eigenvector between neighborhood, the initial value using the class probability value of each pixel on error image as Markov random field model;

S112b, according to maximum a posteriori probability estimation theory and Markov random field theory, construct the markov energy function by described Markov random field model G (V, E) step by step;

According to triumphant Lay Hamiltonian theory, the Markov random field model energy function can be equivalent to characteristic model energy and two part sums of prior model energy.The concrete building method of characteristic model energy and prior model energy is as follows:

1, by the class probability value y of each pixel of Gaussian process sorter output _i ^*As the input feature vector of Markov random field model, suppose that the unique point Output rusults of each node in Markov random field model is labeled as z _iCalculate its characteristic model

Suppose input feature vector

Equal Gaussian distributed, for any one input feature vector y _i ^*, utilize gauss of distribution function, calculate its characteristic model.The gauss of distribution function computing formula that this overdivided region belongs to the k class is as follows:

$p\left(y_i^{*}\right|z_i)=\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\Σ}}_k}}\exp [-\frac{1}{2}{(y_i^{*}-{\mathrm{\μ}}_k)}^T{\mathrm{\Σ}}_k^{-1}(y_i^{*}-{\mathrm{\μ}}_k)]---\left(8\right)$

Wherein, μ _k, ∑ _kThe average and the variance that mean respectively k category feature vector.

2, utilize single order Ising model, calculate the prior model p (z of Markov random field model G (V, E) _i), p (z _i) mean z _iNeighborhood

In pixel between correlativity punishment;

Prior model is taked single order Ising model, only thinks that the mark Output rusults of current certain pixel only is adjacent existence interaction between the element marking Output rusults connect, and and do not have mutual relationship between other element marking Output rusults.Markov single order Ising model formulation is:

$p\left(z_i\right)=\frac{1}{C}\exp (-\underset{j\∈N_{z_i}}{\mathrm{\Σ}}{\mathrm{\β}}_1\·\mathrm{\δ}(z_i,z_j))---\left(9\right)$

Wherein, δ (z _i, z _j) be defined as the Cadillac function, if z _i=z _j, δ (z _i, z _j)=1, otherwise δ (z _i, z _j)=0; β ₁For level and smooth weighting coefficient, control the effect size between the neighborhood overdivided region; C is the Regularization constant.

3, by the characteristic model p (z of each pixel on error image _i| y _i ^*) and prior model p (z _i) summation, the markov energy function of acquisition view picture error image, its function representation form is:

$\mathrm{\Σ}\left(\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\Σ}}_k}}\exp \right[-\frac{1}{2}{(y_i^{*}-{\mathrm{\μ}}_k)}^T{\mathrm{\Σ}}_k^{-1}(y_i^{*}-{\mathrm{\μ}}_k)]+(\frac{1}{C}\exp (-\underset{j\∈N_{z_i}}{\mathrm{\Σ}}{\mathrm{\β}}_1\·\mathrm{\δ}(z_i,z_j)))---\left(10\right)$

Through the simplification to formula (10), its optimized energy function representation form is equivalent to:

$\arg \underset{x\∈X}{\min }\mathrm{\Σ}\{\frac{1}{2}\ln \left(2\mathrm{\π}{\mathrm{\Σ}}_i\right)+\frac{{(y_i^{*}-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{j\∈N_{z_i}}{\mathrm{\Σ}}\mathrm{\δ}(z_i,z_j)\}---\left(11\right)$

S112c step by step, adopt iterated conditional model (iterative conditional mode, be called for short ICM) markov energy function (as shown in Equation 11) is optimized and solves, the general layout of resulting final energy function mark is the result of 2 o'clock phase Remote Sensing Imagery Change Detection.Method is as follows:

1, individual element in error image is carried out to iteration, calculate least energy corresponding to every pixel;

2, in error image the minimum of all pixels can add and, as the energy when time iteration;

3, record the output class mark of the Markov random field model that current least energy is corresponding, be when the suboptimum Output rusults, such mark means that regular rear each pixel of Markov random field model belongs to variation or indeclinable classification logotype;

4,, when the least energy of this iterative computation gained is less than the previous iteration least energy, by general layout corresponding to this associating energy function, i.e. current Markov random field model output class mark set detects Output rusults as changing; Otherwise the output class mark set that last iteration is corresponding detects Output rusults as changing;

5, meet a prior given maximal value when iterations, or, when the iteration Output rusults of continuous 5 times remains unchanged, stop iteration, otherwise, re-execute step 1.In the present invention, the iterations higher limit is set to 200.

Method provided by the invention has proposed a kind of Gaussian process method relevant based on spatial context for Remote Sensing Imagery Change Detection.At first the method utilizes the Gaussian process sorter to obtain initial variation testing result, then utilizes Markov random field model to carry out regular to initial testing result.Characteristics of the present invention are both to have inherited the powerful classification capacity of Gaussian process sorter, can alleviate again the problem that Markov random field relies on initial value, can also utilize the regular item of markov effectively to remove Gauss's sorter and change the noise spot in detecting, improve and change the precision detected.In addition, because the present invention is that to carry out Markov random field model on the initial change testing result basis obtained at the Gaussian process sorter regular, Markov random field model energy function converges faster in actual operation, calculated amount is less, so change and detect with respect to Gaussian process, increase very little computing time.

Below further illustrate application example of the present invention: example has provided one group of typical urban scene QuickBird multispectral image, and the image taking area is the subregion of Beijing, and shooting time is respectively 2002 and 2003, as shown in Figure 3 a and Figure 3 b shows.Image is by (2.4 meters/pixel of full-colour image (0.6 meter/pixel) and multispectral images, red, green, blue, four wave bands of near infrared) form, in order to possess high spatial resolution and spectral resolution simultaneously, full-colour image and multispectral image have been carried out to Ehlers fusion (M.Ehlers.Spectral characteristics preserving image fusion based on Fourier Domain filtering, In Proc.SPIE, 5574:1-4, Bellingham, 2004.), the size of fused image is 1024 * 1024 pixels.Main change in image is: part is waste has built up large stretch of house on the ground, and original construction area has become greenery patches.

The manual detection result of usining as shown in Fig. 4 e, as the reference result, is contrasted method of the present invention and other three kinds of change detecting methods.Wherein, first method is Pixel-level change detecting method based on the Gaussian process sorter (utilize the Gaussian process sorter of describing in the present invention's the 3rd step, the image pixel of take changes detection as unit), and its testing result is as shown in Fig. 4 a.Second method is based on Pixel-level change detecting method (the L. Bruzzone and D.F.Prieto of Markov random field, " Automatic analysis of the difference image for unsupervised change detection; " IEEE Trans.Geosci.Remote Sens.vol.38, no.3, pp.1171-1182, May 2000), its testing result is as shown in Figure 4 b.The third method is based on regional multi-layer change detecting method (F.Bovolo, " A Multilevel Parcel-Based Approach to Change Detection in Very High Resolution " .IEEE Trans.Geosci.Remote Sens.Letters, vol.6, no.1, pp.33-37,2009), its testing result is as shown in Fig. 4 c.The result of utilizing method of the present invention to be detected is as shown in Fig. 4 d.

Estimated (in Table one) from four aspects to changing testing result: 1) false drop rate; 2) loss; 3) error rate; 4) kappa coefficient.

Table one changes testing result and quantizes comparative result

	False drop rate	Loss	Error rate	Kappa
					The first comparative approach	5.27	7.11	12.38	0.68
The second comparative approach	4.44	10.69	15.13	0.59
					The third comparative approach	2.62	8.02	10.65	0.71
Case method of the present invention	7.29	2.55	10.02	0.76

The variation testing result that quantification comparative result explanation the present invention obtains improves a lot than the accuracy of detection of other three kinds of methods.It should be noted that, the above-mentioned definition to each element is not limited in various concrete structures or the shape of mentioning in embodiment, and those of ordinary skill in the art can replace simply to it with knowing, for example:

(1) in the S102 step, watershed segmentation method image over-segmentation method can also be by other image over-segmentation methods, as meanshift method, Ncut method replace, because meanshift method, Ncut method are method well known in the art, be not described in detail herein;

(2) in the S110c step, iterated conditional model (iterative conditional mode is called for short ICM) energy function method for solving can also use figure (Graphcut) method of cutting replace.

Above-described specific embodiment; purpose of the present invention, technical scheme and beneficial effect are further described; institute is understood that; the foregoing is only specific embodiments of the invention; be not limited to the present invention; within the spirit and principles in the present invention all, any modification of making, be equal to replacement, improvement etc., within all should being included in protection scope of the present invention.

Claims (11)

1. 2 o'clock phase method for detecting change of remote sensing image comprise:

Respectively the remote sensing images of 2 o'clock phases are carried out to over-segmentation, obtain a plurality of overdivided regions after dividing processing, the remote sensing images that described 2 o'clock phase remote sensing images are two width the same area different times;

Respectively to per a period of time the phase remote sensing images take pending overdivided region and extract color characteristic, textural characteristics as unit, and entropy feature;

Color, textural characteristics, the entropy feature in corresponding pending zone in 2 o'clock phase remote sensing images are asked to poor, and required difference is as the proper vector of corresponding pending overdivided region in the error image of 2 o'clock phase remote sensing images;

Utilize housebroken Gaussian process sorter, according to the characteristic of correspondence vector, overdivided region in error image is changed class and do not changed two classes classification of class, obtain the class probability of each overdivided region, and by this classification probability assignment to all pixels in this pending overdivided region on error image;

Build Markov random field model on error image, initial value using the class probability value of each pixel in each overdivided region as Markov random field model, calculate the Markov random field model energy function, by energy function optimization, obtain the result of 2 o'clock phase Remote Sensing Imagery Change Detection.

2. 2 o'clock according to claim 1 phase method for detecting change of remote sensing image, describedly utilize housebroken Gaussian process sorter, according to proper vector, overdivided region in error image is changed class and do not change two classes classification of class, and the step of obtaining the class probability of each overdivided region comprises:

Obtain training sample set;

Utilize proper vector and the class mark of the concentrated training sample of training sample to be trained the Gaussian process sorter, obtain the gaussian kernel parameter;

Utilize Gaussian process sorter corresponding to gaussian kernel parameter, pending overdivided region proper vector in described error image is changed class and do not change two classes classification of class, obtain the class probability of each pending overdivided region.

3. 2 o'clock according to claim 1 phase method for detecting change of remote sensing image, the described Markov random field model that builds on error image, initial value using the class probability value of each overdivided region as Markov random field model, calculate the Markov random field model energy function, by energy function optimization, the step of obtaining the result of 2 o'clock phase Remote Sensing Imagery Change Detection comprises:

A Markov random field model G (V of structure on error image, E), in error image, each pixel is expressed as a node V in Markov random field model, the pixel of its peripheral region of each node V forms the neighborhood set N of this node V, the poor E of the limit as Markov random field model of eigenvector between neighborhood, the initial value using the class probability value of each pixel on error image as Markov random field model;

According to maximum a posteriori probability estimation theory and Markov random field theory, by described Markov random field model G (V, E) structure markov energy function;

Adopt the iterated conditional model to be optimized and to solve the markov energy function, the general layout of resulting final energy function mark is the result of 2 o'clock phase Remote Sensing Imagery Change Detection.

4. 2 o'clock according to claim 3 phase method for detecting change of remote sensing image, described according to maximum a posteriori probability estimation theory and Markov random field theory, step by described Markov random field model G (V, E) structure markov energy function comprises:

Input feature vector using the class probability value of each pixel of Gaussian process sorter output as Markov random field model G (V, E), suppose that the unique point Output rusults of each node in Markov random field model is labeled as z _i, calculate its characteristic model

Utilize single order Ising model, calculate the prior model p (z of Markov random field model G (V, E) _i), p (z _i) mean z _iNeighborhood

In pixel between correlativity punishment;

Characteristic model by each pixel on error image

With prior model p (z _i) summation, the markov energy function of acquisition view picture error image.

5. 2 o'clock according to claim 4 phase method for detecting change of remote sensing image, described employing iterated conditional model is optimized and solves the markov energy function, and the step that the general layout of resulting final energy function mark is the result of 2 o'clock phase Remote Sensing Imagery Change Detection comprises:

Step 1, individual element in error image is carried out to iteration, calculate least energy corresponding to every pixel;

The least energy addition of all pixels summation in step 2, error image, as the energy when time iteration;

Step 3, record the output class mark of the Markov random field model that current least energy is corresponding, be when the suboptimum Output rusults, such mark mean Markov random field model regular after each pixel belong to and change or indeclinable classification logotype;

Step 4, when the least energy of this iterative computation gained is less than the previous iteration least energy, by this associating general layout corresponding to energy function, i.e. current Markov random field model output class mark set detects Output rusults as changing; Otherwise the output class mark set that last iteration is corresponding detects Output rusults as changing;

Step 5, the iterations in described step 4 equals M, or, when the inferior iteration Output rusults of N continuous remains unchanged, stops iteration, and the set of current Markov random field model output class mark is detected to Output rusults as changing, otherwise, re-execute step 1.

6. 2 o'clock according to claim 5 phase method for detecting change of remote sensing image, wherein: described M=500, described N=5.

7. according to described 2 o'clock phase method for detecting change of remote sensing image of any one in claim 1 to 6, one of them carries out over-segmentation to the remote sensing images of 2 o'clock phases respectively in the following ways: watershed segmentation method, meanshift method and Ncut method.

8. according to described 2 o'clock phase method for detecting change of remote sensing image of any one in claim 1 to 6, described respectively to per a period of time the phase remote sensing images take pending overdivided region and extract color characteristic, textural characteristics as unit, and also comprise after the step of entropy feature:

2 o'clock each overdivided region corresponding to phase remote sensing images are adjusted, made the overdivided region of correspondence of 2 o'clock phase remote sensing images consistent.

9. according to described 2 o'clock phase method for detecting change of remote sensing image of any one in claim 1 to 6, described 2 o'clock each overdivided region corresponding to phase remote sensing images are adjusted, make in the consistent step of the overdivided region of correspondence of 2 o'clock phase remote sensing images, for an overdivided region A in the first phase image, its corresponding region in second o'clock phase images is overdivided region B, and the adjustment criterion of overdivided region A and overdivided region B is:

If A=B, overdivided region A and overdivided region B all remain unchanged;

If A ≠ B, and A ∩ B=A, overdivided region A remains unchanged, and overdivided region B is divided into to B=B1 ∪ B2, wherein, the corresponding overdivided region that B1 is overdivided region A on second o'clock phase images, B2 is B-B1;

If A ∪ is B>A, and A ∪ B>B, overdivided region A being divided into to A-(A ∩ B) and A ∩ B two parts, overdivided region B is divided into B-(A ∩ B) and A ∩ B two parts.

10. according to described 2 o'clock phase method for detecting change of remote sensing image of any one in claim 1 to 6, described to the time phase remote sensing images take pending overdivided region and extract color characteristic, textural characteristics as unit, and the step of entropy feature comprises:

While take in the phase remote sensing images in the overdivided region target pixel be unit, extract color characteristic, textural characteristics and the entropy feature of each pixel;

Calculate the feature of interior each characteristic mean of pixel of overdivided region as this pending overdivided region.

11. according to described 2 o'clock phase method for detecting change of remote sensing image of any one in claim 1 to 6, described color, textural characteristics, entropy feature by corresponding pending zone in 2 o'clock phase remote sensing images asks poor, and required difference comprises as the step of the proper vector of corresponding pending overdivided region in the error image of 2 o'clock phase remote sensing images:

Respectively color characteristic, textural characteristics, entropy feature are carried out to normalized, obtain normalization color characteristic, textural characteristics and entropy feature afterwards;

According to the over-segmentation image after adjusting, on overdivided region, feature is subtracted each other one by one, and a plurality of features of obtaining pending overdivided region in 2 o'clock phase remote sensing images are poor, the proper vector of corresponding pending overdivided region in the poor formation error image of the plurality of feature.

Images