CN103150718B - Based on the method for detecting change of remote sensing image of Change vector Analysis and classification and predicting - Google Patents

Abstract

Based on a method for detecting change of remote sensing image for Change vector Analysis and classification and predicting, relate to remote sensing images technology, comprising: the first step, respectively color and textural characteristics are extracted in units of pixel to the remote sensing images of two phases; Second step, theoretical with Markov random field, when the two phase Remote Sensing Imagery Change Detection based on classification and predicting method being converted into two, the phase images Er Kefu energy function that starts compares; 3rd step, respectively to twice phase images carry out over-segmentation, respectively two phase over-segmentation image-regions are readjusted according to rule; The two phase Remote Sensing Imagery Change Detection based on Change vector Analysis method are converted to based on provincial characteristics similarity measurement function item by the 4th step; 5th step, combines the markov energy function based on classification and predicting with the similarity measurement function item based on Change vector Analysis, tectonic syntaxis markov energy function; 6th step, solves associating markov energy function with optimization method, exporting change data.

Description

Remote sensing image change detection method based on change vector analysis and classified comparison

Technical Field

The invention relates to the technical field of remote sensing image processing, in particular to a two-time phase remote sensing image change detection method combining a change vector analysis method and a classified comparison method.

Background

Remote sensing image change detection is one of the most important application directions in the field of remote sensing image processing. The remote sensing image change detection technology has important application value in a plurality of military and civil fields such as natural disaster monitoring, national and local resource planning management, military target striking evaluation and the like. With the increase of the resolution of the remote sensing image data, the change information which can be acquired by the change detection is increasingly abundant, so that the actual application range of the change detection is further expanded.

Post-classification comparison (pcc) and change vector analysis (cva) are two commonly used types of change detection methods. The post-classification comparison method generally classifies images of each time phase, and obtains the change of ground coverage in two time phases by directly comparing the labeled images of the classified images. The change vector analysis method directly obtains the change vectors of the two time phase images, and obtains the change situation of the ground coverage on the two time phases through the analysis of the change vectors. The advantages of the post-classification comparison method are: the radiation correction of the images in two time phases is not needed, and the change condition of the coverage type on the ground can be conveniently acquired. The disadvantages are that: neglecting the dependency relationship between the two time phase images easily leads to accumulation of errors, resulting in low detection accuracy. The change vector analysis method has the advantages that: the dependency relationship between the two time phase images is fully utilized, and the accumulation of errors caused by classification errors can be effectively inhibited. The disadvantages are that: detecting the serious influence of radiation factors such as illumination, vegetation and the like; in addition, due to the non-uniqueness of the acquisition of the change vector, the loss of original image information is serious, and the detection result has no uniqueness.

How to comprehensively utilize the advantages of the two methods, realize advantage complementation and improve detection precision is a more meaningful problem. The existing method (J.Chen, X.Chen, X.Cui, and J.Chen, "ChangeVectorrAnalysisinPosterIorProbasilicSpeace: ANewMethodForLandchangeDestection", IEEETrans.Trans.Geosci.RemoteSens.letters, vol.8, No.2, pp.317-321, Mar.2011.) tries to reduce the error accumulation of the classified comparison method by using a change vector analysis method, but the method has the important defects that the dependency information of two-time phase images cannot be fully utilized, manual threshold selection is needed, and the automation degree is not high. The invention comprehensively utilizes the advantages of the comparison method after classification and the change vector analysis method to realize advantage complementation, converts the change detection method combining the comparison method after classification and the change vector analysis method into a combined Markov energy function expression form for the first time, and obtains the optimal change detection output result through the optimization solution of the energy function. The whole detection process does not need manual intervention, and the automation degree is high.

Disclosure of Invention

The invention aims to provide a remote sensing image change detection method based on change vector analysis and classified comparison, so that the advantages of the change vector analysis method and the classified comparison method are fully and effectively utilized, and the change detection precision is improved.

In order to achieve the purpose, the technical scheme of the invention is as follows:

a remote sensing image change detection method based on change vector analysis and classified comparison comprises the following steps:

firstly, extracting color and texture characteristics of the remote sensing images of two time phases by taking a pixel as a unit;

secondly, converting the two-time phase remote sensing image change detection based on the classified comparison method into Markov energy function comparison on the two-time phase image by using a Markov random field theory;

thirdly, respectively over-segmenting the two time phase images, and respectively readjusting the over-segmented image areas of the two time phases according to rules;

fourthly, converting the two-time phase remote sensing image change detection based on the change vector analysis method into a similarity measurement function item based on the regional characteristics;

fifthly, combining the Markov energy function based on classified comparison with a similarity measurement function item based on change vector analysis to construct a combined Markov energy function;

and sixthly, solving the combined Markov energy function by using an optimization method, and outputting a change detection result.

The remote sensing image change detection method comprises the following steps of, in the first step, pixel-level feature extraction, including:

a1, extracting the CIELab color feature, the Gabor texture feature and the entropy feature of each pixel by taking the pixel as a unit; wherein, Lab color characteristic dimension is 3; dimension of the Gabor filter scale parameter and the dimension of the Gabor filter direction parameter are selected according to actual needs, and the dimension of the Gabor texture feature is 5 multiplied by 8 to 40; selecting the size of a window according to the resolution of the image, and calculating the entropy of the image in the window area with the current pixel as the center as the entropy characteristic of the pixel, wherein the dimension of the entropy characteristic is 1;

and a2, respectively carrying out normalization processing on each type of characteristics.

The second step of the remote sensing image change detection method is a post-classification comparison change detection method based on markov energy function comparison: the method comprises the steps of utilizing a Markov random field model to respectively model images of each time phase, constructing a Markov energy function on each time phase image on the basis of the characteristics of the two time phase remote sensing images obtained in the first step, converting the classification problem of the images on each time phase into a Markov energy function optimization problem, and comparing the Markov energy functions of the two time phase images to obtain a two time phase remote sensing image change detection result based on a post-classification comparison method.

The remote sensing image change detection method comprises the following steps:

b1, respectively modeling the image of each time phase according to the maximum posterior probability estimation theory and the Markov random field theory:

a characteristic model: supposing that each type of feature in each time phase image obeys Gaussian distribution, and calculating a feature model by using a Gaussian distribution function; for any pixel p (i, j) in each phase image, the gaussian distribution calculation formula of the pixel belonging to the k-th class is as follows:

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

wherein p (x | y) is a conditional probability form, x represents a feature vector of the current pixel p (i, j), y is a classification output mark corresponding to the current pixel p (i, j), y ∈ { 1.. k }, and mu_k、Σ_kRespectively, mean and variance of the kth class.

A prior model: the prior model adopts a first-order Ising model, namely only the current pixel p (i, j) is considered to have interaction with the pixel of N in the first-order neighborhood, but has no correlation with other pixels; the Markov first-order Ising model is expressed as:

$p\left(y\right)=\frac{1}{Z}\exp (-\underset{y_m\∈N}{\Σ}{\mathrm{\β}}_1\·\mathrm{\δ}\left(y,y_m\right))---\left(2\right)$

wherein p (y) is a prior probability form, (y, y)_m) Defined as the cadilac function: if y is equal to y_mThen (y, y)_m) 1, otherwise (y, y)_m)＝0；β₁Controlling the magnitude of the interaction between the neighborhood pixels for smoothing the weighting coefficients; n is a first-order neighborhood; z is a regularization constant;

b2, Markov energy function construction:

the markov energy on the image of each phase is the sum of the energies of each pixel p (i, j) on the image, and the function expression is as follows:

$\arg \underset{x\∈X}{min}\mathrm{\Σ}\{\frac{1}{2}ln\left(2{\mathrm{\π\Σ}}_k\right)+\frac{{(x-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{y_m\∈N}{\Σ}\mathrm{\δ}(y,y_m)\}---\left(3\right)$

b3, comparing the Markov energy functions of the two-time phase images to obtain a Markov energy function expression form of the classified comparison method: because the images of the two time phases are independent and irrelevant, the classification of each time phase image is independent, and the comparison after the classification of the corresponding two time phase images is that the Markov energy functions of the two time phase images are respectively minimized; directly adding energy functions defined on the two time phase images to be used as an energy function expression form for classifying and comparing the two time phase images:

$\arg \underset{x\∈X}{\min }\underset{s=1}{\overset{2}{\Σ}}\mathrm{\Σ}\{\frac{1}{2}ln\left(2{\mathrm{\π\Σ}}_k\right)+\frac{{(x-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{y_m\∈N}{\Σ}\mathrm{\δ}(y,y_m)\}---\left(4\right)$

the remote sensing image change detection method comprises the third step of adjusting the over-segmentation image area: respectively carrying out over-segmentation processing on the remote sensing images in two time phases, acquiring homogeneous region images with regions as units, and respectively carrying out region adjustment on the two over-segmented images according to rules, wherein the method comprises the following steps:

c1, performing over-segmentation processing on the remote sensing image of each time phase by using a watershed segmentation method, wherein each over-segmentation area is a homogeneous area; when the total number of pixels of a certain area is less than 200, the area stops being divided;

c2, comparing the over-segmentation areas of each time phase image one by one according to the over-segmentation images of the two time phase images, and adjusting the over-segmentation areas on each time phase according to a comparison criterion; assuming that the area a is an over-divided area in the time phase 1 image, and the area B is a corresponding area of the area a on the second time phase image, the comparison criterion between the area a and the area B is as follows:

if A is B, the area A and the area B are kept unchanged;

if A ≠ B and A ≠ B ═ A, then area A remains unchanged, and area B is partitioned into

B-B1 ═ B2, where B1 is the corresponding region of region a on the second phase image and B2 is B-B1;

if A ^ B > A and A ^ B > B, the area A is divided into two parts A- (A ^ B) and A ^ B, and the area B is divided into two parts B- (A ^ B) and A ^ B.

The remote sensing image change detection method, wherein the fourth step is a change detection method based on regional similarity measurement change vector analysis: comparing and analyzing the feature vectors of the two-time phase images, calculating the similarity between the feature vectors, and converting the remote sensing image change vector analysis based on the region as a unit into a similarity measurement function item of the feature vectors, wherein the similarity measurement function item comprises the following steps:

d1, a similarity measure between corresponding pixels p1 and p2 in the two-phase image_pDefined as the characteristic euclidean distance weighted sum of two pixels: the Euclidean distance weighted sum of Lab color features, Gabor texture features and entropy features of pixels p1 and p2, whose function is expressed in the form:

dism_p(p1,p2)＝(x1_Lab-x2_Lab)²+ω₁(x1_Gabor-x2_Gabor)²+ω₂(x1_entory-x2_entory)²(5) wherein, x1_Lab，x1_GaborAnd x1_entoryRespectively representing a Lab color feature vector, a Gabor texture feature and an entropy feature of the pixel p 1; x2_Lab，x2_GaborAnd x2_entoryRespectively representing a Lab color feature vector, a Gabor texture feature and an entropy feature of the pixel p 2; omega₁And ω₂As weighting factors, being constant, typically ω₁∈[0,1]，ω₂∈[0,1]；

d2, similarity measure between any two regions in two-phase images, dis_rI.e. the sum of all pixel similarity measures in the region; the calculation formula is as follows:

dism_r＝∑dism_p(p1,p2)(6)

in the method for detecting a change in a remote sensing image, in the fifth step, a markov energy function is combined to construct: combining the energy items of the change vector analysis with the classified and compared energy functions to obtain a combined Markov energy function; the whole function gives consideration to the fact that the classified comparison method has the optimal solution on one hand and the two time phase images have the maximum similarity on the other hand, so that the change between the two time phase images is optimally detected, and the method comprises the following steps:

adding the energy item based on the change vector analysis method and the energy function based on the classified comparison method to obtain a combined energy expression form based on the change vector analysis method and the classified comparison method:

$\arg \underset{x\∈X}{\min }\underset{s=1}{\overset{2}{\mathrm{\Σ}}}\mathrm{\Σ}\{\frac{1}{2}\ln \left(2{\mathrm{\π\Σ}}_k\right)+\frac{{(x-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{y_m\∈N}{\mathrm{\Σ}}\mathrm{\δ}(y,y_m)+{\mathrm{\β}}_2{\left({\mathrm{dism}}_p\right)}^{\mathrm{\β}}\}---\left(7\right)$

wherein, β₂As weighting factor, the proportion of energy term used for controlling the variation vector analysis method and the energy function of the classified comparison method when β₂When the value is smaller, the comparison method after classification plays a dominant role, when β₂β is a control function, generally in the form of a switch function, and if the scale values corresponding to the pixels in the two time phase images are the same, the value of β is 1, and if the scale values corresponding to the pixels in the two time phase images are different, the value of β is-1.

In the method for detecting a change in a remote sensing image, the sixth step of solving the joint markov energy function by using an optimization method includes:

f1, iterating each pixel of the corresponding coordinates in the two time phase images point by point, and calculating the minimum joint energy of each pixel; summing the minimum joint energy of all pixels in the image as the energy of the current iteration;

f2, recording a change detection output class mark corresponding to the current minimum joint energy, namely the current optimal output result;

f3, iteration: when the minimum joint energy obtained by the current iteration calculation is smaller than the minimum joint energy of the previous iteration, the pattern corresponding to the joint Markov energy function at the current time, namely a current two-time phase image change detection output class diagram is used as an output result; otherwise, the change detection output class diagram corresponding to the last iteration is used as an output result;

f4, when the iteration number satisfies a preset maximum value or when the output result of 5 continuous iterations is kept unchanged, stopping the iteration and outputting the change data.

In the remote sensing image change detection method, the maximum value of the iteration times in the step f4, namely the upper limit value of the iteration times of output change data is set to be 200; the maximum value given in advance by the iteration times is selected in a compromise mode according to the requirements of precision and speed.

In the method for detecting the change of the remote sensing image, iteration is finished in the sixth step, and when the energy function is converged to a local optimal value, the corresponding image pattern is the final change detection result.

The invention provides a change detection method combining change vector analysis and classified comparison aiming at two-time-phase remote sensing images, which respectively converts two-time-phase remote sensing image change detection based on the classified comparison method and two-time-phase remote sensing image change detection based on the change vector analysis method into an energy function expression form by reasonably constructing a Markov random field model energy function, skillfully combines the two energy functions, converts the two-time-phase remote sensing image change detection combining the change vector analysis and the classified comparison into a combined Markov energy function optimization solving problem, and utilizes an iteration condition model to optimize and solve the combined Markov energy function to obtain a change detection result. The method effectively utilizes the advantages of the change vector analysis method and the classified comparison method in the change detection, and improves the precision of the change detection through the advantage complementation, so that the method has better robustness; meanwhile, two-time phase remote sensing image change detection combining change vector analysis and classified comparison is converted into a combined Markov energy function optimization solving problem, so that the automation degree of change detection is improved, and the method has better practicability.

The invention has the advantages of analyzing the change detection method based on the change vector and comparing the change detection method based on the change vector after classification, converts the change detection method combining the comparison method after classification and the change vector analysis method into a combined Markov energy function expression form for the first time, obtains the optimal change detection output result through the optimization solution of the energy function, does not need manual intervention in the detection process and has high automation degree.

Drawings

FIG. 1 is a schematic flow chart of a method for detecting changes in remote sensing images based on analysis and comparison after classification of change vectors according to the present invention;

FIG. 2 is a schematic diagram of a Markov energy function structure for two-time-phase remote sensing image change detection based on a post-classification comparison method in the second step of the method;

FIG. 3 is a schematic diagram of an adjustment criterion of over-segmented image regions in the third step of the method of the present invention;

FIG. 4 is a diagram illustrating the construction of a combined Markov energy function in a fifth step of the method of the present invention;

FIG. 5 shows example data for an application of the method of the present invention; wherein:

fig. 5(a) is a 2002 captured image;

FIG. 5(b) is a 2003 captured image;

FIG. 6 is an exemplary detection result of the method for detecting changes in remote sensing images based on analysis and comparison after classification of change vectors according to the present invention; wherein:

FIG. 6(a) is a variation detection result based on a post-classification comparison method;

FIG. 6(b) is a change detection result based on the change vector analysis method;

FIG. 6(c) is a variation detection result of the combination of variation vector analysis and post-classification comparison in the method of the present invention;

FIG. 6(d) shows the results of the manual detection.

Detailed Description

The invention discloses a remote sensing image change detection method based on change vector analysis and post-classification comparison, which respectively converts two-time phase remote sensing image change detection based on a post-classification comparison method and a change vector analysis method into comparison of Markov energy functions of two-time phase images and similarity measurement energy items, constructs a combined Markov energy function by combining energy expression functions of the two methods, and solves the combined Markov energy function by using the existing energy function optimization solving method to obtain the optimal change detection output result.

Firstly, extracting the basic features of the remote sensing images of two time phases respectively. The method comprises the following steps:

1.1, extracting Lab color characteristics of each pixel by taking the pixel as a unit:

1.1.1 image smoothing filtering: taking the coordinate (i, j) as the center, taking an image area with the size of 3 x 3 (the window size can be adjusted according to the requirement), and calculating the spectral mean value of the pixels in the area as the spectral value of the current pixel p (i, j).

1.1.2 conversion calculation formula from RGB color space to CIELab color space is:

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

$X=\frac{0.412453*R+0.357580*G+0.180423*B}{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)(8)

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.1.3 Lab color feature dimension of the image is 3, and the Lab color feature is stored according to the image coordinate mode.

1.2, extracting Gabor texture features of the image:

1.2.1Gabor Filter:

x'＝xcosθ+ysinθ(9)

y'＝-xsinθ+ycosθ

wherein, λ is wavelength, θ is strip direction of the Gabor filter,the ellipsoid rate of the Gabor is controlled by the phase parameter of the Gabor filter and gamma is the length-width ratio.

The Gabor energy function is defined as follows:

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

wherein,andare respectively g_λ,θ,0(x, y) andthe result of convolution with the grayscale image.

1.2.2 converting the original RGB image into a grayscale image: the spectral values of each pixel in the three channels of image R, G, B are directly averaged together as the spectral value of the grayscale image.

1.2.3 setting the scale parameter of the Gabor filter to be 5 and the direction parameter to be 8, constructing 40 Gabor filters, and filtering the gray image by using each filter, namely convolving each Gabor filter with the gray image; the dimension of the Gabor texture feature is 40;

1.2.4 storing Gabor texture features according to an image coordinate mode;

1.3, extracting entropy characteristics of each pixel by taking pixel coordinates as a unit:

1.3.1 converting the original RGB image into a grayscale image: the spectral values of each pixel in the three channels of image R, G, B are directly averaged together as the spectral value of the grayscale image.

1.3.2 each pixel p (i, j) in the image is centered on the coordinate (i, j), and the entropy of the image area is calculated as the entropy feature of the pixel, taking the image area of 10 × 10. The formula for calculating entropy is as follows:

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

wherein p is_iWhich represents the proportion of the number of pixels with a gray value i in the image area to the number of pixels in the entire image area.

1.3.3 dimension of entropy feature of image is 1, storing entropy feature according to image coordinate mode.

And 1.4, normalizing each type of feature respectively. The normalized calculation formula is:

f＝(f-min(min(f)))/(max(max(f))-min(min(f)))(12)

and secondly, respectively modeling the images of each time phase by using a Markov random field model, constructing a Markov energy function on each time phase image respectively on the basis of the obtained remote sensing image feature vectors of the two time phases, comparing the Markov energy functions of the two time phases, and obtaining a Markov energy function expression form which is compared after classification. The method comprises the following steps:

2.1, constructing a Markov random field model G (V, E) on the image of each time phase based on the images of the two time phases: each pixel p (i, j) in each time phase image is represented as a node V in the Markov random field model, each node V and four pixels with the Euclidean distance difference value of 1 around the node V form a first-order neighborhood N of the node V, and the characteristic vector difference between the pixels in 4 vertical and horizontal directions in the neighborhood is used as an edge E of the Markov random field model.

And 2.2, classifying the remote sensing images of each time phase respectively on the basis of the feature vectors extracted from the images of each time phase. The classification principle is as follows: pixels with higher similarity in the image are grouped into one class so as to ensure that the intra-class variance is smaller than the inter-class variance, and the smaller the ratio is, the better the ratio is.

2.2.1 according to the maximum posterior probability estimation theory, giving a feature vector X of a certain time phase image, and obtaining a maximum posterior probability solving problem by the classification output Y:

p(Y|X)＝p(X|Y)p(Y)/p(X)(13)

2.2.2 corresponding to each pixel p (i, j) on the image, the mathematical description of the classification process is expressed as:

$y=\arg \underset{x\∈X}{max}\left\{p\left(x\right|y)p\left(y\right)\right\}---\left(14\right)$

wherein, x is the characteristic vector of the current pixel p (i, j), and y is the classification mark corresponding to the point; p (x | y) describes the distribution relation between the image features and the classification output class labels; and p (y) describing the interaction relationship between the output class labels of the images.

2.3 according to the Markov random field model-maximum posterior probability estimation theory, constructing a characteristic model p (x | y) and a prior model p (y) on the image of each time phase:

2.3.1 feature model: and (4) calculating a feature model by using a Gaussian distribution function on the assumption that each type of feature in each phase image obeys Gaussian distribution. For any pixel p (i, j) in each phase image, the gaussian distribution function calculation formula of the pixel belonging to the k-th class is as follows:

$p\left(x\right|y)=\frac{1}{\sqrt{2{\mathrm{\π\Σ}}_k}}\exp \[-\frac{1}{2}{(x-{\mathrm{\μ}}_k)}^T{\mathrm{\Σ}}_k^{-1}(x-{\mathrm{\μ}}_k)\]---\left(15\right)$

wherein x represents the feature variation vector of the current pixel p (i, j), y is the classification output mark corresponding to the current pixel p (i, j), y ∈ { 1.. k }, μ_k、Σ_kRespectively representing the mean and variance of the kth class of feature vectors.

2.3.2 prior model: the prior model takes a first-order Ising model, that is, only the current pixel p (i, j) is considered to have interaction with the pixels of N in the first-order neighborhood, but has no interaction relation with other pixels. The Markov first-order Ising model is expressed as:

$p\left(y\right)=\frac{1}{Z}\exp (-\underset{y_m\∈N}{\Σ}{\mathrm{\β}}_1\·\mathrm{\δ}\left(y,y_m\right))---\left(16\right)$

wherein, (y, y)_m) Defined as a cadilac function if y ═ y_mThen (y, y)_m) 1, otherwise (y, y)_m)＝0；β₁Controlling the magnitude of the interaction between neighboring pixels for smoothing the weighting coefficients(ii) a N is a first-order neighborhood; z is a normalization constant.

2.4 constructing a Markov energy function on each phase image respectively: the markov energy in the image at each time phase is the sum of the energies of each pixel p (i, j) in the image. And (3) respectively taking logs at two sides of the posterior probability expression, so that the maximum posterior probability equivalence is expressed as the following functions:

$\arg \underset{x\∈X}{min}\mathrm{\Σ}\{\frac{1}{2}ln\left(2{\mathrm{\π\Σ}}_k\right)+\frac{{(x-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{y_m\∈N}{\Σ}\mathrm{\δ}(y,y_m)\}---\left(17\right)$

when the class mark value of each point in the image is most suitable, the value of the whole energy function tends to be stable and has the minimum value, and the corresponding Markov energy pattern is the optimal classification output result at the moment.

2.5, since the images of the two phases are independent and uncorrelated, the classification of each phase image is independent and uncorrelated. When the Markov energy function on the image of each phase is minimum, the optimal output class diagram of the image corresponding to each phase is obtained. Thus, a post-classification comparison corresponding to two-phase images is equivalent to a minimization of the sum of the Markov energy functions of the two-phase images. The energy functions defined on the two time phase images are directly added to be used as an energy function expression form for classifying and comparing the two time phase images:

$\arg \underset{x\∈X}{min}\underset{s=1}{\overset{2}{\Σ}}\mathrm{\Σ}\{\frac{1}{2}ln\left(2{\mathrm{\π\Σ}}_k\right)+\frac{{(x-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{y_m\∈N}{\Σ}\mathrm{\δ}(y,y_m)\}---\left(18\right)$

and thirdly, respectively carrying out over-segmentation processing on the remote sensing images of the two time phases to obtain a homogeneous region image with a region as a unit, carrying out region comparison on the two over-segmented images according to a rule, and respectively readjusting the regions of the two time-phase over-segmented images. The method comprises the following steps:

3.1, by means of a watershed segmentation method (J.roerdinkanda.Meijster, "the Watershed transforms: Definitions, Algorithms, and ParallellizationStrategies," FundamentaInformationae, vol.41, pp.187-228,2000.), the remote sensing images of each phase are respectively subjected to segmentation processing based on the CIELab color of each phase image. Each over-segmentation region has certain region homogeneity consistency, namely a homogeneous region. In the over-segmentation process, in order to ensure the maximum consistency of the region and the actual interpretability of change detection, the area of the region is selected as a constraint condition, namely when the total number of pixels of a certain segmented region is less than 200, the region stops segmentation;

and 3.2, comparing the over-segmentation areas of each time phase image one by one according to the over-segmentation images of the two time phase images, and adjusting the over-segmentation areas on each time phase according to a comparison criterion. Assuming that the area a is an over-divided area in the time phase 1 image, and the area B is a corresponding area of the area a on the second time phase image, the comparison criterion between the area a and the area B is as follows:

if A is B, the area A and the area B are kept unchanged;

if A ≠ B and A ≠ B ═ A, then area A remains unchanged, and area B is partitioned into

B-B1 ═ B2, where B1 is the corresponding region of region a on the second phase image and B2 is B-B1;

if A ^ B > A and A ^ B > B, the area A is divided into two parts A- (A ^ B) and A ^ B, the area B is divided into two parts B- (A ^ B) and A ^ B;

and fourthly, comparing and analyzing the feature vectors of the two time phase images, calculating the similarity between the feature vectors, and converting the remote sensing image change vector analysis based on the region as a unit into the similarity measurement of the feature vectors. The method comprises the following steps:

4.1 similarity metric function dis of pixels p1 and p2_p: similarity measure dism between arbitrary pixels p1 and p2 in two-phase images_pDefined as the euclidean distance weighted sum of the two pixel features, whose function is expressed as:

dism_p(p1,p2)＝(x1_Lab-x2_Lab)²+ω₁(x1_Gabor-x2_Gabor)²+ω₂(x1_entory-x2_entory)²，(19)

wherein, x1_Lab，x1_GaborAnd x1_entoryRespectively representing a Lab color feature vector, a Gabor texture feature and an entropy feature of the pixel p 1; x2_Lab，x2_GaborAnd x2_entoryRespectively representing a Lab color feature vector, a Gabor texture feature and an entropy feature of the pixel p 2; omega₁And ω₂As weighting factors, being constant, typically ω₁∈[0,1]，ω₂∈[0,1]. Usually, omega is selected according to the characteristics of the image₁And ω₂The value of (c). In the invention, omega₁And ω₂The values are all 1.

4.2 similarity measure between any two regions in two-phase images, dism_rI.e. the sum of all pixel similarity measures in the region. The calculation formula is as follows:

dism_r＝∑dism_p(x1,x2)(20)

and fifthly, combining the energy item for analyzing the change vector with the classified and compared energy function to obtain a combined Markov energy function. The whole function ensures that the classified comparison method has the optimal solution on one hand and the similarity of the two time phase images is maximum on the other hand, so that the change between the two time phase images is optimally detected. The specific method is that a similarity measurement function item of a change vector analysis method is used as a biased item of an energy function of a classified comparison method, and a new energy function expression is constructed:

5.1, adding the energy item of the change vector analysis method and the energy function after classification comparison to obtain a combined Markov energy function expression form of the change vector analysis method and the classified comparison method:

$\arg \underset{x\∈X}{\min }\underset{s=1}{\overset{2}{\Σ}}\mathrm{\Σ}\{\frac{1}{2}ln\left(2{\mathrm{\π\Σ}}_k\right)+\frac{{(x-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+{\mathrm{\β}}_1\underset{y_m\∈N}{\Σ}\mathrm{\δ}(y,y_m)+{\mathrm{\β}}_2{\left({\mathrm{dism}}_p\right)}^{\mathrm{\β}}\}---\left(21\right)$

wherein, β₂As weighting factor, the proportion of energy term used for controlling the variation vector analysis method and the energy function of the classified comparison method when β₂When the value is smaller, the comparison method after classification plays a dominant role, when β₂β is a control function, generally in the form of a switch function, and if the scale values corresponding to the pixels in the two time phase images are the same, the value of β is 1, and if the scale values corresponding to the pixels in the two time phase images are different, the value of β is-1.

5.2 in the present invention, β₁Take the constant 1, β₂Take the constant 0.5.

And sixthly, solving the combined Markov energy function by using an optimization method. The method comprises the following steps:

6.1, solving the energy function (21): and (2) adopting an Iterative Condition Model (ICM) to carry out optimization solution on the combined Markov energy function (21). The number of iterations is typically around 200, eventually converging to a local optimum.

6.1.1 assume that in the post-classification comparison method, the number of classes of image classification for each phase is k, and the value of k is generally given directly.

6.1.2 respectively calculating Markov energy when the class mark value of each time phase image is [1,2, …, k ], namely calculating the value of a time function (17) when the class mark value is [1,2, …, k ] for each pixel in the two time phase images;

6.1.3 calculating a similarity measure between corresponding pixels of the two time phase images, namely calculating the value of a function term (19);

6.1.4 selecting the value of the switching function beta according to the class mark values of the two time phase pixels; adding Markov energy of a classified comparison method of pixels in the two-time phase image and Markov energy of a change vector analysis method to obtain combined Markov energy on the current pixel, namely energy of each pixel in a function (21);

6.1.5 calculating the joint Markov energy of all pixels in the image, summing to obtain the joint Markov energy of the two-time phase image, namely a function value corresponding to the function (21);

6.1.6 recording the class label graph of the two time phase image corresponding to the output energy value of the current function (21);

6.1.7 compare the two phases of the object-like map pixel by pixel: if the pixel class values corresponding to the two time phase images are the same, the change detection output mark corresponding to the pixel position is 0; otherwise, the pixel position output is marked as 1;

6.1.8 iteration: when the value of the combined Markov energy function (21) obtained by the current iteration calculation is smaller than the previous iteration value, the pattern corresponding to the current combined Markov energy function, namely the current two-time phase image change detection output class diagram is used as an output result; otherwise, the change detection output class diagram corresponding to the last iteration is used as an output result;

6.1.9 the iteration is stopped when the number of iterations satisfies a predetermined maximum value or when the output of 5 consecutive iterations remains unchanged. In the method, the upper limit value of the iteration times is set to be 200 (the iteration times can be selected according to the compromise between the precision and the speed requirement).

6.1.10, the iteration is finished, and the corresponding image pattern is the final change detection result when the energy function converges to the local optimum value.

Examples of applications of the invention are further illustrated below:

the example shows a set of typical urban scene Quickbird multispectral images, the image shooting region is a partial region of Beijing, and the shooting time is 2002 and 2003 respectively. The image consists of a panchromatic image (0.6 m/pixel) and a multispectral image (2.4 m/pixel, four bands of red, green, blue and near infrared), and in order to have high spatial resolution and spectral resolution simultaneously, the panchromatic image and the multispectral image are subjected to Ehlers fusion (m.ehlers. spectral characteristics prediction imaging fusion based on the image analysis of the blue and multispectral images, inproc. spie,5574:1-4, Bellingham,2004.), and the size of the fused image is 1024 × 1024 pixels. The main changes in the image are: a large number of houses are built on a part of wasteland, and the original building area becomes a green land.

And comparing the accuracy of the three change detection methods by taking the manual detection result as a reference result, wherein one change detection method is a change detection method based on a comparison method after classification, one change detection method is a change detection method based on a change vector analysis method, and the other change detection method is the embodiment method of the invention. Wherein, the change detection method based on the classified comparison method is the output result of the second step in the specific implementation steps of the invention; the change detection method based on the change vector analysis method is to directly analyze the change vector by using a Markov random field (L.BruzzoneandD.F.Prieto, "Automatic natriented information for unsupervised change detection," IEEETrans. Geosci. RemoteSens. vol.38, No.3, pp.1171-1182, May2000.)

The results of the change detection were evaluated in four ways (see table 1): 1) the false detection rate; 2) the omission rate; 3) an error rate; 4) the kappa coefficient.

TABLE 1 Change detection results quantitative comparison results

The quantitative comparison result shows that the detection precision of the remote sensing image change detection method based on change vector analysis and classified comparison is greatly improved.

Claims (8)

1. A remote sensing image change detection method based on change vector analysis and classified comparison is characterized by comprising the following steps:

firstly, extracting color and texture characteristics of the remote sensing images of two time phases by taking a pixel as a unit;

secondly, converting the two-time phase remote sensing image change detection based on the classified comparison method into Markov energy function comparison on the two-time phase image by using a Markov random field theory;

thirdly, respectively over-segmenting the two time phase images, and respectively readjusting the over-segmented image areas of the two time phases according to rules;

fourthly, converting the two-time phase remote sensing image change detection based on the change vector analysis method into a similarity measurement function item based on the regional characteristics;

fifthly, combining the Markov energy function based on classified comparison with a similarity measurement function item based on change vector analysis to construct a combined Markov energy function;

solving the combined Markov energy function by using an optimization method, and outputting a change detection result;

wherein, in the second step, the detection method of the comparison change after classification based on the Markov energy function comparison comprises the following steps: respectively modeling images of each time phase by using a Markov random field model, respectively constructing a Markov energy function on each time phase image on the basis of the characteristics of the two time phase remote sensing images acquired in the first step, converting the classification problem of the images on each time phase into a Markov energy function optimization problem, and comparing the Markov energy functions of the two time phase images to acquire a change detection result of the two time phase remote sensing images based on a post-classification comparison method;

in the third step, the two time phase images are respectively over-segmented, and the readjustment of the two time phase over-segmented image areas according to the rule comprises the following steps:

c1, performing over-segmentation processing on the remote sensing image of each time phase by using a watershed segmentation method, wherein each over-segmentation area is a homogeneous area; when the total number of pixels of a certain area is less than 200, the area stops being divided;

c2, comparing the over-segmentation areas of each time phase image one by one according to the over-segmentation images of the two time phase images, and adjusting the over-segmentation areas on each time phase according to a comparison criterion; assuming that the area a is an over-divided area in the time phase 1 image, and the area B is a corresponding area of the area a on the second time phase image, the comparison criterion between the area a and the area B is as follows:

if A is B, the area A and the area B are kept unchanged;

if a ≠ B, and a ≠ B ═ a, then area a remains unchanged, partitioning area B into B ═ B1 ═ B2, where B1 is the corresponding area of area a on the second phase image and B2 is B-B1;

if A ^ B > A and A ^ B > B, the area A is divided into two parts A- (A ^ B) and A ^ B, and the area B is divided into two parts B- (A ^ B) and A ^ B.

2. The remote sensing image change detection method according to claim 1, wherein the first step, pixel-level feature extraction, comprises:

a1, extracting Lab color features, Gabor texture features and entropy features in a CIELab color space of each pixel by taking the pixel as a unit; wherein, Lab color characteristic dimension is 3; dimension of the Gabor filter scale parameter and the dimension of the Gabor filter direction parameter are selected according to actual needs, and the dimension of the Gabor texture feature is 5 multiplied by 8 to 40; selecting the size of a window according to the resolution of the image, and calculating the entropy of the image in the window area with the current pixel as the center as the entropy characteristic of the pixel, wherein the dimension of the entropy characteristic is 1;

and a2, respectively carrying out normalization processing on each type of characteristics.

3. The method for detecting changes in remote sensing images according to claim 1, wherein the second step specifically comprises:

b1, respectively modeling the image of each time phase according to the maximum posterior probability estimation theory and the Markov random field theory:

a characteristic model: supposing that each type of feature in each time phase image obeys Gaussian distribution, and calculating a feature model by using a Gaussian distribution function; for any pixel p (i, j) in each phase image, the gaussian distribution calculation formula of the pixel belonging to the k-th class is as follows:

wherein p (x | y) is a conditional probability form, and x represents a feature vector of the current pixel p (i, j); y is the current imageThe classification output label, y ∈ { 1.. multidot.k }; mu, corresponding to element p (i, j)_k、Σ_kRespectively representing the mean and the variance of the kth class;

a prior model: the prior model adopts a first-order Ising model, namely only the current pixel p (i, j) is considered to have interaction with the pixel of N in the first-order neighborhood, but has no correlation with other pixels; the Markov first-order Ising model is expressed as:

wherein p (y) is a prior probability form, (y, y)_m) Defined as the cadilac function: if y is equal to y_mThen (y, y)_m) 1, otherwise (y, y)_m)＝0；β₁Controlling the magnitude of the interaction between the neighborhood pixels for smoothing the weighting coefficients; n is a first-order neighborhood; z is a regularization constant;

b2, Markov energy function construction:

the markov energy on the image of each phase is the sum of the energies of each pixel p (i, j) on the image, and the function expression is as follows:

wherein X is a characteristic vector of the time phase image;

b3, comparing the Markov energy functions of the two-time phase images to obtain a Markov energy function expression form of the classified comparison method: because the images of the two time phases are independent and irrelevant, the classification of each time phase image is independent, and the comparison after the classification of the corresponding two time phase images is that the Markov energy functions of the two time phase images are respectively minimized; directly adding energy functions defined on the two time phase images to be used as an energy function expression form for classifying and comparing the two time phase images:

4. a remote sensing image change detection method according to claim 3, characterized in that in the fourth step, the change detection method based on the analysis of the regional similarity measure change vector: comparing and analyzing the feature vectors of the two-time phase images, calculating the similarity between the feature vectors, and converting the remote sensing image change vector analysis based on the region as a unit into a similarity measurement function item of the feature vectors, wherein the similarity measurement function item comprises the following steps:

d1, a similarity measure between corresponding pixels p1 and p2 in the two-phase image_pDefined as the characteristic euclidean distance weighted sum of two pixels: the Euclidean distance weighted sum of Lab color features, Gabor texture features and entropy features of pixels p1 and p2, whose function is expressed in the form:

dism_p(p1,p2)＝(x1_Lab-x2_Lab)²+ω₁(x1_Gabor-x2_Gabor)²+ω₂(x1_entory-x2_entory)²，(5)

wherein, x1_Lab，x1_GaborAnd x1_entoryRespectively representing a Lab color feature vector, a Gabor texture feature and an entropy feature of the pixel p 1; x2_Lab，x2_GaborAnd x2_entoryRespectively representing a Lab color feature vector, a Gabor texture feature and an entropy feature of the pixel p 2; omega₁And ω₂As a weighting coefficient, being a constant, ω₁∈[0,1]，ω₂∈[0,1]；

d2, similarity measure between any two regions in two-phase images, dis_rI.e. the sum of all pixel similarity measures in the region; the calculation formula is as follows:

dism_r＝∑dism_p(p1,p2)(6)。

5. the remote sensing image change detection method according to claim 4, wherein in the fifth step, a joint Markov energy function is constructed: combining the energy items of the change vector analysis with the classified and compared energy functions to obtain a combined Markov energy function; the whole function gives consideration to the fact that the classified comparison method has the optimal solution on one hand and the two time phase images have the maximum similarity on the other hand, so that the change between the two time phase images is optimally detected, and the method comprises the following steps:

adding an energy item based on a change vector analysis method and an energy function based on a classified comparison method to obtain a combined Markov energy expression form based on the change vector analysis method and the classified comparison method:

wherein, β₂As weighting factor, the proportion of energy term used for controlling the variation vector analysis method and the energy function of the classified comparison method when β₂When the value is smaller, the comparison method after classification plays a dominant role, when β₂β is a control function in the form of a switch function, if the scale values corresponding to the pixels in the two time phase images are the same, the value of β is 1, and if the scale values corresponding to the pixels in the two time phase images are different, the value of β is-1.

6. A remote sensing image change detection method as claimed in claim 1, wherein said sixth step of solving the joint markov energy function using an optimization method comprises:

f1, iterating each pixel of the corresponding coordinates in the two time phase images point by point, and calculating the minimum joint energy of each pixel; summing the minimum joint energy of all pixels in the image as the energy of the current iteration;

f2, recording a change detection output class mark corresponding to the current minimum joint energy, namely the current optimal output result;

f3, iteration: when the minimum joint energy obtained by the current iteration calculation is smaller than the minimum joint energy of the previous iteration, the pattern corresponding to the joint Markov energy function at the current time, namely a current two-time phase image change detection output class diagram is used as an output result; otherwise, the change detection output class diagram corresponding to the last iteration is used as an output result;

f4, when the iteration number satisfies a preset maximum value or when the output result of 5 continuous iterations is kept unchanged, stopping the iteration and outputting the change data.

7. A remote sensing image change detection method as claimed in claim 6, characterized in that the maximum value of the iteration number in the step f4, namely the upper limit value of the iteration number of the output change data, is set to 200; the maximum value given in advance by the iteration times is selected in a compromise mode according to the requirements of precision and speed.

8. A remote sensing image change detection method as recited in claim 6, characterized in that, when iteration is finished in the sixth step and the energy function converges to the local optimum value, the corresponding image pattern is the final change detection result.