CN100538399C - Cutting apart and classification and variety detection integration method of high-resolution remote sensing image - Google Patents

Abstract

The invention discloses a kind of high-resolution remote sensing image that develops based on multilevel collection cuts apart and classification and variety detection integration method.It is characterized in that: (1) image pre-service (radiation, registration and characteristic filtering); (2) set up cutting apart and disaggregated model that multilevel collection develops, determine the initial profile of each level set function automatically, the 1st o'clock phase image cut apart and classified with GIS data behind the registration; (3) still adopt model described in (2), optimize the initial profile of each level set function, the phase image all adopts increment type to cut apart and classify to the 2nd to T the time; (4) be the unit to cut apart the back object, image classification result mutually when relatively two of i and i+1 are adjacent determines region of variation; (5) be back to (3), the phase image cuts apart and classification and change-detection when finishing all T.Advantage: compare towards the K averaging method of pixel with tradition, classification and accuracy of detection all increase, and are suitable for the change-detection of sequence remote sensing image, at aspects such as disaster monitoring and survey of territorial resources wide application are arranged.

Description

Cutting apart and classification and variety detection integration method of high-resolution remote sensing image

Technical field

The invention belongs to computing machine remote sensing image information extraction technology field, relate to a kind of high spatial resolution remote sense image that develops based on multilevel collection curved surface and cut apart and classification and variety detection integration method.

Background technology

The remote sensing image information extraction is a basic and important research contents in fields such as computing machine and remote sensing.But along with sensor technology, satellite communication technology, rapid development of computer technology, remote sensing image more and more presents the characteristics of " more than three ", i.e. multisensor, multiple space or spectral resolution, multidate, and this makes remote sensing image automated information retrieval difficulty more.

In recent years, (spatial resolution is less than or equal to 5 meters to high spatial resolution, hereinafter to be referred as high resolving power) satellite, as IKONOS, QuickBird, and domestic military satellite, moonlet etc. all have very high spatial resolution, and it provides abundant terrestrial object information for the user, but has brought bigger difficulty also for the automated information retrieval task.Along with high resolution image is increasingly extensive in the application of aspects such as land resources survey, city planning, disaster monitoring, traditional no longer suitable towards the cutting apart of pixel, classification and change detecting method, the deficiency of its performance mainly is: noise and too much details usually can cause cutting apart, " spiced salt " phenomenon occur among classification and the change-detection result in (1) image, its output result can't represent a significant target, does not match with actual atural object; (2) it does not take into full account the texture of high resolution image except that spectrum, the features such as spatial relationship between structure, geometry and atural object, and nicety of grading is low; (3) it is cut apart with classification results and is not easy to direct vector quantization, is not easy to set up with geographic information data get in touch.Previous work has proved that OO information extraction technology is to alleviate these not enough effective means, pertinent literature has: G..J.Hay, G.Castilla, " Object-based image analysis:strengths; Weakness; Opportunities and Threats (SWOT); " OBIA2006:International Archive of Photogrammetry, Remote Sensing and Spatial Information Sciences, July, 2006.http: //www.commission4.isprs.org/obia06/papers.htm.Multi-scale division among the commercial packages eCognition and in the classification representative of cutting apart on the basis of this class technology just, it improved that tradition is cut apart and classification results in " spiced salt " phenomenon, improved nicety of grading.In addition, the classification on watershed divide (watershed), the average drifting multi-scale division bases such as (meanshift) also is an effective way of handling high-resolution remote sensing image.But in these class methods, the performance of dividing method has a significant impact nicety of grading.

Change-detection is an important content in remote sensing image information extraction technology field, a plurality of links such as it relates to image radiant correction, registration, cuts apart, Classification and Identification, and error propagation is serious, and present change detecting method roughly divides two classes:

● image algebraic operation method.It is very high to radiant correction and registration accuracy requirement, is not suitable for the change-detection problem of multidate Heterogeneous Sensor image.

● classification back relative method.In this method, classification is crucial, again it is divided into towards pixel and OO change detecting method according to the sorting technique difference, the former nicety of grading is low, classification results is not easy to the GIS database update, and very high to the registration error requirement, is generally less than 0.5 pixel; Latter's nicety of grading height, registration error is required relatively low, error amount seeing image image space resolution and change-detection target and decide.In addition, these class methods are mostly at the remote sensing image of limited number of phases when a plurality of (2-10 time phase).

For the high-resolution remote sensing image, adopt OO thought to carry out that image is cut apart or classification or change-detection are a kind of effective way, referring to document: G.Willhauck, Comparison of object oriented classificationtechniques and standard image analysis for the use of change detection between SPOTmultispectral satellite images and aerial photos.ISPRS, Vol.XXXIII, Amsterdam, 2000, and Niemeyer, I.﹠amp; M.J.Canty, Pixel-Based and Object-Oriented Change DetectionAnalysis Using High-Resolution Imagery.Proc.25th Symposium on Safeguards andNuclear Material Managment, Stockholm, Sweden, 13-15 May 2003.Based on above thought, this patent discloses a kind of new OO high-resolution remote sensing image information extracting method, and it is a kind of based on cutting apart and classification and variety detection integration method that multilevel collection curved surface develops.

Level set theory (level set theory) is to handle sealing moving interface effective computational tool of how much change in topology in the evolutionary process in time.Its ultimate principle is: by level set function Φ (x, y, t) Biao Shi three-dimension curved surface is in image, under the effect of external force, change its topology, when curved surface evolution speed goes to zero, the plane curve of the implicit expression of the zero level collection of this level set function can be represented target or zone boundary, pertinent literature has: S.Osher, J.A Sethian, " Fronts Propagatingwith Curvature Dependent Speed:Algorithms Based on Hamilton-Jacobi Formulations, " Journalof Computational Physics, vol.79, pp.12-49,1988 and S.Osher, R.Fedkiw.Level Set Methods andDynamic Implicit Surfaces.New York:Springer-Verlag New York Inc., ch.3, pp.25-37,2003.It changes into three-dimensional curved surface evolution problem with the image segmentation problem of two dimension, by finding the solution the image segmentation that realizes of the partial differential equation (group) relevant with time t.It is a kind of driving wheel contour method that topology freely changes that has, and can realize the image division under the global energy optiaml ciriterion.Its advantage is: range statistics feature or the goal gradient feature of not only considering image, and utilized the geometrical property (curvature, normal force) of curve, and can merge priori of target etc., at last these information are integrated in the EVOLUTION EQUATION of level set function as curve evolvement power.

Since two thousand, the level set theory becomes the focus of field scholar's research such as state, inside and outside fluid mechanics, materialogy, mathematics, computer vision, Flame Image Process gradually.Wherein, the research in medical image and use a lot ofly, pertinent literature has: L.Vese, T.F.Chan, " A Multiphase Level Set Framework for Image Segmentation Using theMumford and Shah Model, " International Journal of Computer Vision, vol.50, no.3, pp.271-293, July, 2002, and document H.K Zhao, T.F.Chan, B.Merriman, S.Osher, " Variational Level SetApproach to Multiphase Motion; " Journal of Computational Physics, no.0167, pp.179-195, February, 1996.This has proved the feasibility of level set theory in image segmentation and classification, but has its deficiency in the above method, is subjected to the restriction of level set number or has the fuzzy phenomenon etc. of cutting apart as number of categories.

Be different from medical image, remote sensing image has characteristics such as diversity, complexity of imaging environment and multiband, multidate, the high spatial resolution of imaging platform, and therefore, propose effective information extracting method is the focus that people pay close attention to always.Multilevel collection theoretical remote sensing image cut apart or aspect research such as classification seldom.Pertinent literature has: A.R.Mansouri, A.Mitiche, C.Va ' zquez, " Multiregion competition:A level set extension of regioncompetition to multiple region image partitioning; " Computer Vision and ImageUnderstanding, no.101, pp.137-150, in 2006, the author mainly is at medical science, the dividing method that the multilevel collection that video image proposes develops, but do not relate to high spatial resolution or multispectral image, and do not take into full account the imaging circumstances complexity of remote sensing image, characteristics such as type of ground objects is various.Document: O.Besbesy, Z.Belhadj, N.Boujema, " AdaptiveSatellite Images Segmentation by Level Set Multiregion Competition; " INRIA, Technical Report, no.5855, March, 2006 work based on A.R.Mansouri have proposed an adaptive disaggregated model; Document: Ball, J.E.Bruce, L.M.Level set segmentation of remotely sensed hyperspectral images.IGARSS 2005.vol.8, page (s): among the 5638-5642, the author has set up the disaggregated model of multilevel collection at Hyperspectral imaging, but it underuses the texture of remote sensing image, several how more characters of ground object, accuracy still remains to be improved.Method in the above pertinent literature all is supervision property, needs sample learning and parameter estimation before cutting apart or classifying, and has increased workload;

Application facet in change-detection, people such as the Vivien Mallet of Florida State University have developed the Multivac software package based on the thought that single level set develops, be used for the fire spreading process is carried out the growth of modeling and simulation silicon nano material, pertinent literature is seen http://vivienmallet.net/fronts/index.php.In addition, document Dell ' Acqua, F.Gamba, PPrevedini, P., Level-set based extraction and tracking of meteorological objects in satellite images, IGARSS 2000.vol.2, page (s): among the 627-629, utilize the dynamic evolution properties of level set theory to design the cloud Target Tracking System of middle low resolution weather satellite data, more than work and confirmed level set theory applicability aspect the single goal change-detection in many multi-temporal remote sensing images, but utilize the variation monitoring of the multi-class targets that relates in the tasks such as status investigation at the soil, above method is invalid, needs to adopt multilevel collection (the level set function number is greater than 1) to cut apart or sorting technique.

In Flame Image Process, cut apart similar with the notion of classification.In fields such as computer visions, generally claim image segmentation.Because existing level set is cut apart the researcher and is mainly come from applied mathematics, computer vision and field of medical image processing, therefore in the pertinent literature in this field, usually use this term of image segmentation, do not relate to the notion of classification; And in pattern-recognition and remote sensing field, cutting apart with classification is two different processes, conceptive also have a difference, the former is meant image is divided into a plurality of homogeneous regions, be indifferent to the number in homogeneity district and the atural object attribute in each homogeneity district, and the latter is divided into different classes of atural object with image, and the general user is concerned about the classification number that divided and the atural object attribute of each class very much.Level set develop theoretical remote sensing image cut apart and classification application research in, part author represents assorting process with cutting apart term, as the document that author A.R.Mansouri et.al. and author O.Besbesy et.al. deliver, in fact, they have only realized the process of image classification.The present invention is towards the remote sensing field, it has distinguished this two notions, the method that is proposed has not only realized classification, and provided image segmentation result, but cut apart different with the expression of classification results, the former be with the symbol of each level set function with whole image be divided into different classes of after, demarcate all enclosed region that comprise of all categories more one by one, form the homogeneity district; Be divided into whole image different classes of and the latter only is a symbol with each level set function.

Summary of the invention

Problem to be solved by this invention is: provide a kind of high-resolution remote sensing image that develops based on multilevel collection curved surface to cut apart and the integral method of classification and change-detection, the evolution that this method is object with the zone by the closed curve family in plane realizes cutting apart of image and classifies, be a kind of new OO cutting apart and classification and change-detection thought, effectively alleviated traditional " spiced salt " phenomenon in the information extracting method of pixel.Compare with the sorted change detecting method of traditional K average, its classification and change-detection precision are higher.

Technical scheme provided by the invention is: a kind of high-resolution remote sensing image that develops based on multilevel collection curved surface is cut apart and the integral method of classification and change-detection, may further comprise the steps:

One, the pre-service of high-resolution remote sensing image:

The phase high-resolution remote sensing image carries out radiant correction respectively during a) to each; And the phase high-resolution remote sensing image carries out registration during to behind the radiant correction each, and then with the GIS data registration of areal;

B), form m dimensional feature vector image to high-resolution remote sensing image texture feature extraction and spectral signature after radiant correction and the registration

Wherein, m is the dimension of the eigenvector image of extraction, and m 〉=1; u _kThe characteristic image of representing k component, 1≤k≤m; Utilize the multiple dimensioned incorgruous diffusion technique of vector total variation to carry out filtering then;

Two, foundation is cut apart and disaggregated model to the eigenvector image after the pre-service:

At first, set up multilevel collection evolution equations, i.e. partial differential equations:

$\left\{\begin{array}{c}\frac{{\∂\mathrm{\Φ}}_1}{\∂t}(x,y)=[-{\mathrm{\ξ}}_1(x,y)+{\mathrm{\φ}}_1(x,y)+{\mathrm{\μCur}}_{{\mathrm{\Φ}}_1}]\left|\right|{\▿\mathrm{\Φ}}_1(x,y)\left|\right|\\ \·\\ \·\\ \·\\ \frac{{\∂\mathrm{\Φ}}_i}{\∂t}(x,y)=[-{\mathrm{\ξ}}_i(x,y)+{\mathrm{\φ}}_i(x,y)+{\mathrm{\μCur}}_{{\mathrm{\Φ}}_i}]\left|\right|{\▿\mathrm{\Φ}}_i(x,y)\left|\right|\\ \·\\ \·\\ \·\\ \frac{{\∂\mathrm{\Φ}}_{N-1}}{\∂t}(x,y)=[-{\mathrm{\ξ}}_{N-1}(x,y)+{\mathrm{\φ}}_{N-1}(x,y)+{\mathrm{\μCur}}_{{\mathrm{\Φ}}_{N-1}1}]\left|\right|{\▿\mathrm{\Φ}}_{N-1}(x,y)\left|\right|\end{array}\right.---\left(1\right)$

In the formula (1),

Be implicit expression N-1 dimension level set function vector, N is the atural object classification number that high-resolution remote sensing image comprises; Function ξ _i(x y) is i classification C among the Ω of bidimensional image territory _iStatistical nature, computing formula is:

${\mathrm{\ξ}}_i(x,y)=\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{(u_k(x,y)-c_{\mathrm{ik}})}^4,1\≤i\≤N---\left(2\right)$

In the formula (2), c _IkBe i classification C on the k dimensional feature vector image among the Ω of bidimensional image territory _iAverage, computing formula is:

$c_{\mathrm{ik}}=\frac{{\∫}_{C_i}{u}_k(x,y)\mathrm{dxdy}}{{\∫}_{C_i}\mathrm{dxdy}}---\left(3\right)$

In the formula (3), integral sign

The expression limit of integration is i classification C _iThe pixel set that comprises; In the formula (1), function phi _i(x, computing formula y) is:

${\mathrm{\φ}}_i(x,y)={\mathrm{\ξ}}_{i+1}(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{i+1}(x,y,t)>0\}}(x,y)+{\mathrm{\ξ}}_{i+2}(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{i+1}(x,y,t)\≤0\}}(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{i+2}(x,y,t)>0\}}(x,y)+\·\·\·$

$+{\mathrm{\ξ}}_{N-1}(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{i+1}(x,y,t)\≤0\}}(x,y)\·\·\·{\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{N-2}(x,y,t)\≤0\}}(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{N-1}(x,y,t)>0\}}(x,y)---\left(4\right)$

$+{\mathrm{\ξ}}_N(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{i+1}(x,y,t)\≤0\}}(x,y)\·\·\·{\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{N-2}(x,y,t)\≤0\}}(x,y){\mathrm{\χ}}_{\{{\mathrm{\Φ}}_{N-1}(x,y,t)\≤0\}}(x,y)$

In the formula (1), Be the mean curvature of i level set function, computing formula is: ${\mathrm{Cur}}_{{\mathrm{\Φ}}_i}=\▿\·\frac{{\▿\mathrm{\Φ}}_i}{\left|\right|{\▿\mathrm{\Φ}}_i\left|\right|},$

Be gradient operator; Parameter μ is the weights of mean curvature;

Be the indicative function of i level set function, be defined as:

Specifying constraint:

$C_i(x,y,t)=\underset{1\≤i\≤N-1}{\min }\left\{{\mathrm{\ξ}}_i(x,y,t)\right\}---\left(5\right)$

In the formula (5), C _i(x, y, t) in the t time iteration of expression, ξ _i(x, y, t) pixel (x, class label y) of i level set function correspondence of value minimum; When multilevel collection evolutionary process satisfies the condition of convergence, class label C _iVariation remain static;

The partial differential equations (1) that has constraint condition (5) is exactly cutting apart of high-resolution remote sensing image and disaggregated model, and separate (being that image is cut apart and classification results) of this model represented in the following ways: the 1st is expressed as to the N-1 class: C _i(x, y, t)=(x, y) ∈ Ω | Φ _i(x, y, t)〉0}, and wherein, 1≤i≤N-1, the N class is expressed as: $C_N(x,y,t)=\{(x,y)\∈\mathrm{\Ω}|\underset{i=1}{\overset{N-1}{\mathrm{}\∩\mathrm{}}}{\mathrm{\Φ}}_i(x,y,t)\≤0\},$ And the zone that segmentation result is exactly the zero level collection of all level set functions to surround;

Three, the cutting apart and classify of phase high-resolution remote sensing image the 1st time:

A) utilize method of finite difference to carry out discretize to cutting apart with disaggregated model, boundary condition is the Neumann boundary value condition; Parameter in input type (1)-(3): the atural object classification number N that image comprises gets positive integer, and N 〉=2, and like this, the required level set function number of cutting apart with disaggregated model is N-1; The weights μ of mean curvature=(0.1) ^p* 255 ², p is an integer, span is 0～5; Spatial sampling interval delta x=1 on the row, column direction, Δ y=1, i.e. a pixel; Iteration time step delta t depends on the evolution speed of spatial sampling interval and level set function, and its value must satisfy CFL condition;

B) in the high-resolution remote sensing image level set of the 1st o'clock phase being cut apart and is classified, cut apart with the starting condition method to set up of disaggregated model as follows: the GIS of this area data behind registration are down auxiliary, the initial profile of each level set function is set, concrete grammar: utilize the locus of each atural object classification in the GIS data and attribute to come initial circle of picture, and with each the initial round starting condition of symbolic distance conversion of carrying out as each level set function, adopt in the step 3 a) described boundary condition and discretize mode and parameter value realize the 1st o'clock mutually the level set of image cut apart and classify, cut apart and obtain cutting apart when classify end after plane closure curve family;

Four, the high-resolution remote sensing image of phase adopts increment type to cut apart and classification the 2nd to T the time: mutually high-resolution remote sensing image still adopts cutting apart and disaggregated model described in the step 2 to the 2nd to T the time, but adopt new starting condition and new image action scope, new starting condition method to set up is as follows: the closed curve family in plane after the phase high-resolution remote sensing image is cut apart during with j is the phase high-resolution remote sensing image preceding initial profile of cutting apart and classify during as the j+1 that is adjacent; Still adopt in the step 3 a) when described boundary condition and discretize mode and parameter value are realized j+1 cutting apart and classifying of mutually high-resolution remote sensing image, the closed curve family in plane after obtaining cutting apart when finishing cut apart and classify; Wherein: T be the change-detection image the time mutually the sum; J is an integer variable, the high-resolution remote sensing image of phase when representing j, and span is for being 1≤j≤T-1;

Five, the change-detection of high-resolution remote sensing image mutually when two of j and j+1 are adjacent: from j=1, when j during with j+1 mutually high-resolution remote sensing image cut apart finish with assorting process after, begin to these two that the phase high-resolution remote sensing image carries out change-detection when adjacent, concrete grammar is: the classification results that contrasts these two high-resolution remote sensing images of phase when adjacent, with wherein any one the time phase high-resolution remote sensing image classification results be reference, another the time phase high-resolution remote sensing image classification results in, with the object after cutting apart is comparing unit, if change has taken place in the whole pixels of any one object or the class label of part pixel, variation has taken place in the type of ground objects that then shows this object, and the locus and the category attribute of output all changes pixel;

Six, returned for the 4th step, the process of circulation execution in step four to six, the phase image cuts apart and classification and change-detection task when having finished all T.

To size be more than 512 * 512 T not simultaneously the high-resolution remote sensing image of phase cut apart and classification and change-detection can adopt following steps: one carries out pre-service set by step; Adopt the wavelet multiresolution rate to decompose again; Cut apart and classification and change-detection according to above-mentioned steps two to six, up to having finished all T cutting apart and classification and change-detection task of phase image simultaneously not.

Under the situation of the GIS polar plot that does not have areal, the starting condition method to set up of cutting apart with disaggregated model described in the above-mentioned steps two is: draw the initial circle of any radius automatically equably in image action scope scope, carry out the symbolic distance conversion again and come each level set function of initialization.

The present invention has set up a kind of multilevel collection evolutionary model based on the zone and has realized that OO image is cut apart and classification and change detecting method, and its advantage shows:

A) but it utilizes the multi-source information amalgamation of level set framework, has set up fusion spectrum, texture, several how manifold multiwave high-resolution remote sensing image and has cut apart and disaggregated model; It is the simplification to A.R.Mansouri proposition method, parameter is few and need not sample learning and parameter estimation procedure, having optimized and cut apart and classification speed, is the QuickBird image of 512 * 512 four wave bands to size, and it is cut apart and expend time in average out to 77 seconds of assorting process;

B) realized OOly cutting apart and classifying, avoided " spiced salt " phenomenon towards the pixel method, alleviated traditional watershed divide, over-segmentation that shows in the region growing method or less divided phenomenon, nicety of grading has improved 2%-10%, compare traditional sorted change detecting method of K average, its accuracy of detection improves 1%-8%, and the information extraction result is easy for the renewal of GIS database;

C) it utilizes free topology change characteristics such as division that the level set curved surface shows, merging in evolutionary process, realized cutting apart and classification of high-resolution remote sensing image simultaneously, improved and cut apart the pattern of independently carrying out in the existing method with classification, reduce the error propagation approach, do not needed the aftertreatment of classifying;

D) it relies on continuity and the dynamic that the level set curved surface develops, the thought that has adopted increment type to cut apart and classify in the classification link of change-detection; It quickens the cutting apart and assorting process of high-resolution remote sensing image of phase when single by the initialization condition of Optimization Model, compare with the change-detection process of not optimizing starting condition, this increment type is cut apart with sorting technique makes the whole spent time of change-detection process shorten 5%-15%;

E) the present invention is the OO back change detecting method of cutting apart and classify, and registration error is had certain robustness.

The present invention can effectively reduce empty inspection and false drop rate in military target detects; Utilize aspect the change-detection in disaster surveillance, soil, this method is flexible, and model parameter is less; Compare traditional sorted change detecting method of K average, to the accuracy of detection raising 1%-8% of sequence remote sensing image, and the information extraction result is easy for the geographic information database renewal.

Compare with the multiresolution dividing method in existing K average and the eCognition software package, the deficiency that still has the speed aspect in the information extraction of this method to the significantly remote sensing image multiclass atural object more than 512 * 512, but aspect the numerical solution of equation, also have the potentiality of optimizing, therefore very big application prospect is arranged in sensor information extractive technique field.

Description of drawings

Fig. 1 is a master routine operational flow diagram of the present invention.

Fig. 2 is cut apart and classify end time domain and classification demarcation synoptic diagram of the present invention.

Fig. 3 is the increment type of the present invention synoptic diagram of cutting apart and classify.

Fig. 4 is the OO change-detection synoptic diagram in classification of the present invention back.

Embodiment

Referring to Fig. 1, Fig. 2, Fig. 3 and Fig. 4, the invention provides a kind of based on cutting apart and classification and variety detection integration method that multilevel collection develops.If the time number of phases of image is T (T for greater than 1 positive integer), then to the 1st to T the time phase cut apart and classification and change-detection implementing procedure referring to accompanying drawing 1, comprise the steps:

1, the pre-service of raw video

A) common, absolute radiation is proofreaied and correct desired parameters and is difficult to obtain, therefore, to the homogeneity sensor not simultaneously the phase remote sensing image adopt common relative radiant correction, as histogram matching method, dark collection-bright collection method etc., these two kinds of methods reach more multi-method referring to document: Ding Lixia, Zhou Bin, royal people's tide; 5 kinds of relative radiation correction method researchs in the remote sensing monitoring. journal of Zhejiang university (agricultural and life science version) .2005,31 (3): 269～276; To the not phase remote sensing image simultaneously of Heterogeneous Sensor, under the situation that does not possess the absolute radiation correction parameter, omit the radiant correction operation;

B) high-resolution remote sensing image behind radiant correction utilization polynomial method or rational function method etc. are carried out interworking standard between the image, concrete operation method is referring to document: Sun Jiabing, oxazepan, Guan Zequn, remote sensing principle, methods and applications, Beijing: Mapping Press, 1997; And document: Zhang Yongsheng, the superfine work of Gong Dan, high-definition remote sensing satellite application---imaging model, Processing Algorithm and application technology, Beijing: Mapping Press, 2005;

C) to the not image mutually simultaneously after radiant correction and the registration, [concrete operation method is referring to document: Zhao Yindi to adopt Gabor wave filter or nonlinear organization tensor method texture feature extraction, the expansion of Markov model and the application in the high resolution image classification thereof, the doctorate paper, Wuhan University, 2006, and Thomas Brox, Joachim Weickerta.A TV flow based local scale estimate and its application to texture discrimination, Journalof Visual Communication and Image Representation, 2006,17 (5): 1053-1073], (the feature selecting concrete operation method is referring to document: Sun Jiabing for spectral band after all of comprehensive raw video or the feature selecting again, oxazepan, Guan Zequn, remote sensing principle, methods and applications, Beijing: Mapping Press, 1997), together constitute m dimensional feature vector space, wherein, m is the dimension of the eigenvector of extraction, and m 〉=1; Utilize vector total variation stream to carry out feature space filtering again, improve the heterogeneity of intra-zone, strengthen interregional feature difference, and keeping important zone boundary, concrete operation method is referring to document: P.Blomgren, T.F.Chan, " Color TV:Total VariationMethods for Restoration of Vector-Valued Images; " IEEE Transactions on ImageProcessing, 1998,7 (3): 304-309.

2, setting up the image that multilevel collection develops cuts apart and disaggregated model

A) foundation is about the energy equation of the closed curve family in plane

Supposing to be positioned on the bidimensional image action scope Ω one is the closed curve bunch of parameter with arc length s, time t ${\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_i(x(s,t),y(s,t)),1<i<N-1,$ Set up curve family ${\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_i(x(s,t),y(s,t)):\left[\mathrm{0,1}\right]\&RightArrow;\mathrm{\&Omega;},$ If zone

Represent curve respectively

Inside and outside zone, the energy equation of curve family can be expressed as formula (1), the concrete connotation of this equation and method for building up are referring to document: A.R.Mansouri, A.Mitiche, C.Va ' zquez, " Multiregioncompetition:A level set extension of region competition to multiple region imagepartitioning; " Computer Vision and Image Understanding, no.101, pp.137-150, October, 2006.

$E\left[{\left\{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_i\right\}}_{i=1}^{N-1}\right]={\&Integral;}_{R_{{\stackrel{\&RightArrow;}{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}}_1}}{\mathrm{\&xi;}}_1(x,y)\mathrm{dxdy}+{\&Integral;}_{R_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_2}}{\mathrm{\&xi;}}_2(x,y)\mathrm{dxdy}+...$

$+{\&Integral;}_{R_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_1}^c\&cap;R_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_2}^c\&cap;...\&cap;R_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{k-1}}^c\&cap;R_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_k}}{\mathrm{\&xi;}}_k(x,y)\mathrm{dxdy}+...$

$+{\&Integral;}_{R_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_1}^c\&cap;R_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_2}^c\&cap;...\&cap;R_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{N-2}}^c\&cap;R_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{N-1}}^c}{\mathrm{\&xi;}}_N(x,y)\mathrm{dxdy}$ (1)

In the formula (1), N is the atural object classification number that predefined image comprises; Last is the penalty term of curve family total length, and μ is a penalty factor, and μ is constant and μ〉0, this penalty term has guaranteed the slickness of partitioning boundary, and this penalty factor also can be considered cuts apart scale parameter, and this value is big more, and the cutting object area is big more.It is worth μ=(0.1) ^p* 255 ², p is an integer, span is 0～5;

Energy equation (1) has utilized N-1 to tie up parameterized curve family

Image territory Ω is divided into N classification, that is, $\mathrm{\&Omega;}=\underset{i=1}{\overset{N}{\mathrm{}\&cup;\mathrm{}}}{C}_i,$ C _iRepresent i classification, each classification all is a pixel set.Because high resolution image is cut apart and is one with assorting process and continues to optimize the process of finding the solution, therefore, C _iPixel number that comprises and locus all constantly change in time, when energy reaches overall minimum or part hour, class label C _iVariation remain static, cut apart with assorting process finish this moment.Although this energy equation (1) form is identical with energy equation in the described method of author A.R.Mansouri, its concrete calculation expression is improved.Improvements mainly are:

I) ξ _i(x, y) connotation of function and calculating aspect: in the method, ξ _i(x, y) function has reflected the zone

Statistical nature, it is the main drive of entire curve bunch evolution, shown in its computing formula is shown suc as formula (2).Energy theorem (1) after improving is primarily aimed at cutting apart of multiband high-resolution remote sensing image and classifies, as the image after the panchromatic and multispectral fusion of IKONOS; It has not only considered the spectral signature of high-resolution remote sensing image, has also made full use of its important spatial texture feature, make cut apart with classification results more accurate; In addition, ξ _i(x, calculating y) need not class probability distribution hypothesis, thereby has saved loaded down with trivial details sample learning and parameter estimation procedure, makes energy equation (1) parameter still less, calculates easier; Importantly, ξ _i(x, y) these computing method of function have been avoided covariance matrix and inversion process thereof, and having adopted power is 4 data degree of confidence item, has accelerated the evolution of curve family.And.Attention: in the following elaboration, the energy equation (1) after energy equation (1) is all represented to improve among the present invention.

${\mathrm{\&xi;}}_i(x,y)=\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{(u_k(x,y)-c_{\mathrm{ik}})}^4,1\&le;i\&le;N---\left(2\right)$

In the formula (2),

Be filtered eigenvector image, m is the dimension of the eigenvector image of extraction, u _kThe characteristic image of representing k component, 1≤k≤m; c _IkBe i classification C on the k dimensional feature vector image in the bidimensional image territory _iThe average of all pixels that comprised is along with the increase of iterations, c _IkValue is tending towards the average of actual area, and its computing formula (3) is:

$c_{\mathrm{ik}}=\frac{{\&Integral;}_{C_i}{u}_k(x,y)\mathrm{dxdy}}{{\&Integral;}_{C_i}\mathrm{dxdy}}---\left(3\right)$

In the formula (3), integral sign

The expression limit of integration is i classification C _iThe pixel set that comprises;

Ii) the present invention has provided the constraint condition that blur prevention is cut apart, and its expression formula is suc as formula (4):

$C_i(x,y,t)=\underset{1\&le;i\&le;N-1}{\min }\left\{{\mathrm{\&xi;}}_i(x,y,t)\right\}---\left(4\right)$

The implication of formula (4) is: in each iteration step, select ξ _i(x, y), the class label of i the level set function correspondence that 1≤i≤N-1 value is minimum is as t moment pixel (x, y) final class label C _iWhen multilevel collection evolutionary process satisfies the condition of convergence, class label C _iVariation remain static.This constraint condition eliminated this method to high-resolution remote sensing image cut apart with classification results in the fuzzy phenomenon of cutting apart.

B) energy equation (1) utilization first variation method is derived EVOLUTION EQUATION about curve family.

Because the energy equation (1) that provides in the 2nd step in the embodiment is identical with the energy equation that author A.R.Mansouri proposes in form, therefore the derivation of curve family EVOLUTION EQUATION has all utilized the first variation principle to minimize the energy functional of each curve, utilizes the gradient descent method to draw the EVOLUTION EQUATION of curve family again.Specific operation process is referring to document: A.R.Mansouri, A.Mitiche, C.Va ' zquez, " Multiregion competition:A level set extension ofregion competition to multiple region image partitioning; " Computer Vision and ImageUnderstanding, no.101, pp.137-150,2006 and document: A.R.Mansouri, A.Mitiche, C.Va ' zquez, " Multiregion Competition:A Level Set Extension of Region Competition to MultipleRegion Partitioning of Images and Image Sequences; " Computer Vision and ImageUnderstanding, 2004.Like this, draw curve family

EVOLUTION EQUATION, i.e. partial differential equations shown in the formula (6):

$\left\{\begin{array}{c}\frac{\&PartialD;{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_1}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_1\left({\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_1(x,y)\right)+{\mathrm{\&phi;}}_1\left({\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_1(x,y)\right)+\mathrm{\&mu;}{\mathrm{Cur}}_1(x,y)]{\stackrel{\&RightArrow;}{n}}_1(x,y)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \frac{\&PartialD;{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_i}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_i\left({\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_i(x,y)\right)+{\mathrm{\&phi;}}_i\left({\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_i(x,y)\right)+\mathrm{\&mu;}{\mathrm{Cur}}_i(x,y)]{\stackrel{\&RightArrow;}{n}}_i(x,y)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \frac{\&PartialD;{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{N-1}}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_{N-1}\left({\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{N-1}(x,y)\right)+{\mathrm{\&phi;}}_{N-1}\left({\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{N-1}(x,y)\right)]+\mathrm{\&mu;}{\mathrm{Cur}}_{N-1}(x,y)]{\stackrel{\&RightArrow;}{n}}_{N-1}(x,y)\end{array}\right.---\left(5\right)$

In the formula (5),

Cur _iRepresent curve respectively

Unit outside method vector and mean curvature.And

${\mathrm{\&phi;}}_i(x,y)={\mathrm{\&xi;}}_{i+1}(x,y){\mathrm{\&chi;}}_{{\stackrel{\&RightArrow;}{R}}_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{i+1}}}(x,y)+{\mathrm{\&xi;}}_{i+2}(x,y){\mathrm{\&chi;}}_{{\stackrel{\&RightArrow;}{R}}_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{i+1}^c}}(x,y){\mathrm{\&chi;}}_{{\stackrel{\&RightArrow;}{R}}_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{i+2}}}(x,y)+\&CenterDot;\&CenterDot;\&CenterDot;$

$+{\mathrm{\&xi;}}_{N-1}(x,y){\mathrm{\&chi;}}_{{\stackrel{\&RightArrow;}{R}}_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{i+1}^c}}(x,y)\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&chi;}}_{{\stackrel{\&RightArrow;}{R}}_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{N-2}^c}}(x,y){\mathrm{\&chi;}}_{{\stackrel{\&RightArrow;}{R}}_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{N-1}}}(x,y)$

$+{\mathrm{\&xi;}}_N(x,y){\mathrm{\&chi;}}_{{\stackrel{\&RightArrow;}{R}}_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{i+1}^c}}(x,y)\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&chi;}}_{{\stackrel{\&RightArrow;}{R}}_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{N-2}^c}}(x,y){\mathrm{\&chi;}}_{{\stackrel{\&RightArrow;}{R}}_{{\stackrel{\&RightArrow;}{\mathrm{\&gamma;}}}_{N-1}^c}}(x,y)$

C) derivation of level set EVOLUTION EQUATION

The level set theory is found the solution problem to the curve shown in equation (5) and had the following advantages: the curve topology need not parametrization when freely changing, and the zero level collection of plane closed curve available horizontal collection curved surface is represented.If the closed curve family in plane

Corresponding level set function vector is

The closed curve family in plane is corresponding one by one with each component of level set function vector.Therefore, plane closed curve

The available horizontal set function

The zero level collection ${\stackrel{\&RightArrow;}{\mathrm{\&Phi;}}}_i=0$ Expression, and

Inside and outside zone

With

Correspond respectively to the zone ${\stackrel{\&RightArrow;}{\mathrm{\&Phi;}}}_i>0$ With ${\stackrel{\&RightArrow;}{\mathrm{\&Phi;}}}_i<0,$ Then equation (5) is transformed into the partial differential equations (6) that level set function is represented

$\left\{\begin{array}{c}\frac{{\&PartialD;\mathrm{\&Phi;}}_1}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_1(x,y)+{\mathrm{\&phi;}}_1(x,y)+{\mathrm{\&mu;Cur}}_{{\mathrm{\&Phi;}}_1}]\left|\right|{\&dtri;\mathrm{\&Phi;}}_1(x,y)\left|\right|\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \frac{{\&PartialD;\mathrm{\&Phi;}}_i}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_i(x,y)+{\mathrm{\&phi;}}_i(x,y)+{\mathrm{\&mu;Cur}}_{{\mathrm{\&Phi;}}_i}]\left|\right|{\&dtri;\mathrm{\&Phi;}}_i(x,y)\left|\right|\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \frac{{\&PartialD;\mathrm{\&Phi;}}_{N-1}}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_{N-1}(x,y)+{\mathrm{\&phi;}}_{N-1}(x,y)+{\mathrm{\&mu;Cur}}_{{\mathrm{\&Phi;}}_{N-1}1}]\left|\right|{\&dtri;\mathrm{\&Phi;}}_{N-1}(x,y)\left|\right|\end{array}\right.---\left(6\right)$

In the formula (6): function phi _i(x, y) expression formula is:

${\mathrm{\&phi;}}_i(x,y)={\mathrm{\&xi;}}_{i+1}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)>0\}}(x,y)+{\mathrm{\&xi;}}_{i+2}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)\&le;0\}}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+2}(x,y,t)>0\}}(x,y)+\&CenterDot;\&CenterDot;\&CenterDot;$

$+{\mathrm{\&xi;}}_{N-1}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)\&le;0\}}(x,y)\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-2}(x,y,t)\&le;0\}}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-1}(x,y,t)>0\}}(x,y)$

$+{\mathrm{\&xi;}}_N(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)\&le;0\}}(x,y)\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-2}(x,y,t)\&le;0\}}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-1}(x,y,t)\&le;0\}}(x,y)$

In the following formula,

The mean curvature of representing i level set function, its expression formula is: ${\mathrm{Cur}}_{{\mathrm{\&Phi;}}_i}=\&dtri;\&CenterDot;\frac{{\&dtri;\mathrm{\&Phi;}}_i}{\left|\right|{\&dtri;\mathrm{\&Phi;}}_i\left|\right|}$

Be gradient operator); Be the indicative function of i level set function, be defined as:

The partial differential equations (6) that has constraint condition (4) is exactly cutting apart of high-resolution remote sensing image and disaggregated model.Separate (being that image is cut apart and classification results) of this model represented in the following ways: the 1st is expressed as to the N-1 class: C _i(x, y, t)=(x, y) ∈ Ω | Φ _i(x, y, t)〉0}, and wherein, 1≤i≤N-1, the N class is expressed as: $C_N(x,y,t)=\{(x,y)\&Element;\mathrm{\&Omega;}|\underset{i=1}{\overset{N-1}{\mathrm{}\&cap;\mathrm{}}}{\mathrm{\&Phi;}}_i(x,y,t)\&le;0\},$ And the zone that segmentation result is exactly the zero level collection of all level set functions to surround;

This method for expressing is different from the described method of author A.R.Mansouri.In addition, the described method of A.R.Mansouri is mainly used in computer vision field, and the author, is not cut apart but relate to so author's method is that remote sensing image has been implemented classification in fact with cutting apart the notion that image classification in the remote sensing field expressed in this term.

3, the 1st o'clock phase image cuts apart and classifies

Image cut apart and classification is exactly a process of finding the solution the numerical solution of the partial differential equations (6) that has constraint condition (4).Formula (6) be one have just, the partial differential equations under the continuous meaning of boundary value condition.In the discrete two dimensional image, adopt numerical solution to approach truly and separate, its process that quantizes is as follows:

A) adopt this simple and rapid discretize of method of finite difference mode, wherein, the space local derviation is to approach with central difference method, and the time local derviation adopts the forward direction eulerian difference method to approach; Spatial sampling interval delta x=1 on the row, column direction, Δ y=1, i.e. a pixel; Iteration time step delta t depends on the evolution speed of spatial sampling interval and level set function, and its value must satisfy CFL condition, to keep the stability of evolutionary process; Adopt the heavy initial method of level set function of second order accuracy, guarantee higher accuracy and moderate calculation cost.Level set function is defined as a symbolic distance function, and as the starting condition of partial differential equations (6), purpose is to keep the slickness of level set curved surface and the convenience that the zero level collection is found the solution, its initial method as the formula (7):

In the formula (7), (x y) is the Euclidean distance function to d; Boundary condition is the Neumann boundary value condition, and purpose is to make the level set curved surface at image boundary

The place stops to develop, and its expression formula is formula (8):

$\frac{{\&PartialD;\mathrm{\&Phi;}}_i(x,y,t)}{\&PartialD;{\stackrel{\&RightArrow;}{n}}_i}=0$ When $(x,y,t)\&Element;\&PartialD;\mathrm{\&Omega;}\&times;[0,+\&infin;]---\left(8\right)$

Formula (8) expression: the normal derivative at image boundary place level set function is zero.About the Euclidean distance function, CFL condition constraint value and the level set function of Δ t down weighs the specific operation process of initial method referring to document: S.Osher, R.Fedkiw.Level Set Methods and Dynamic Implicit Surfaces.New York:Springer-VerlagNew York Inc., ch.3, pp.25-37,2003;

B) setting of the starting condition of partial differential equations (6) is very important, and its locus in image and topology have very big influence to the evolution speed of curve family.Specifically, final cutting apart and classification results more approached in the position of initial profile, and then required time cost is few more; And initial profile can be the plane closed curve of arbitrary shape in principle, but its topology more approaches target to be split or region shape, and then curve family convergence required time is short more, promptly cuts apart with the assorting process consumed time few more.Traditional level set cut apart with sorting technique in all adopt manual type to set the initial profile of level set function, automaticity is low.The GIS data, as topomap, thematic map, it has comprised physical message, shape size, spatial relationship of ground object target etc., the distribution status that has reflected atural object largely, it can assist remote sensing image to cut apart and classify, and more detailed content is referring to document: Wang Ying, the Liu Min warbler, Huang Wenqian; GIS is to the booster action of classification of remote-sensing images interpretation; Marine charting, 2002,22 (3): 12-15.Therefore, being provided with of high-resolution remote sensing image starting condition for the 1st o'clock phase utilized the GIS data of the areal of registration, its priori that is considered as classification is instructed the distribution of initial profile, improve automaticity and the accuracy of cutting apart and classifying with this.Concrete grammar is: utilize the locus of existing atural object classification in the GIS data and attribute to come initial circle of picture, and with each the initial round starting condition of symbolic distance conversion of carrying out as each level set function, the concrete grammar of symbolic distance conversion is referring to document S.Osher, R.Fedkiw.Level Set Methods and DynamicImplicit Surfaces.New York:Springer-Verlag New York Inc., ch.3, pp.25-37,2003.Do not having under the areal GIS data conditions, the high-resolution remote sensing image of first o'clock phase is cut apart and the starting condition method to set up of classifying is: the initial circle of drawing any radius in image action scope scope automatically equably, carry out the symbolic distance conversion again and come each level set function of initialization, like this, automaticity has improved, but this method is responsive to the locus of initial profile, cuts apart with nicety of grading to be lower than precision when not adopting the auxiliary initialization of GIS data.

C) condition of convergence of cutting apart and classify: cut apart with assorting process be exactly at the beginning of, under boundary value condition (7) and the constraint condition (4), find the solution iteratively and artificial relevant partial differential equations (6) of time, when the closed curve family convergence in plane, cut apart with assorting process and finish, the condition of convergence and the parameter that relates to are referring to document: H.K Zhao, T.F.Chan, B.Merriman, S.Osher, " Variational Level Set Approach to Multiphase Motion, " Journal of ComputationalPhysics, no.0167, pp.179-195, February, 1996.

D) cut apart expression with classification results: the separating of partial differential equations (6) that has constraint condition [formula (4)] is one group of curve family that is positioned at imagery zone and class border; Adopt the 2nd step of concrete implementation content c) in the expression mode of model solution of elaboration to representing the cutting apart of high-resolution remote sensing image of the 1st o'clock phase with classification results.But the prerequisite that adopts this expression mode is a partial differential equations (6) must possess the described constraint condition of formula (4), and the effect of this constraint condition is to avoid the fuzzy phenomenon of cutting apart and classify, and avoids promptly that (x y) locates, Φ at pixel _i(x, y, t)〉0 and Φ _s(x, y, t)〉0 (1≤i, s≤N-1, i ≠ s).The implication of this expression mode is: when curve family is restrained, if pixel (x, y) symbol of locating i level set function is for just, and then the class label of this pixel is C _i, but the N class is by not belonging to the 1st set that the pixel of any class constitutes to the N-1 class among the Ω of image territory.And the zone that final segmentation result is exactly the curve at the zero level collection place of all level set functions to surround.See shown in the accompanying drawing 2 that three level set functions are divided into four classifications with image territory Ω, i.e. N=4, cut apart with classification results in, the class label is C ₁Comprising a class label in the top left region at place is C ₂Homogeneity district A ₁, the class label is C ₂Comprising a class label in the right regions at place is C ₃Homogeneity district A ₂, the class label is C ₃Comprising a class label in the zone, lower-left at place is C ₄Homogeneity district A ₃, C ₄The lower right area at place does not comprise any homogeneity district.Be similar to A ₁, A ₂, A ₃Such homogeneity district is exactly the object after cutting apart, and when N got other value, the segmentation result recognition methods was identical therewith.More than set forth and shown that multilevel collection evolutionary model can obtain to cut apart and classification results simultaneously.At last, model according to step 2 foundation, adopt initial and boundary value condition and the discretize mode described in the step 3 to carry out numerical solution, realize cutting apart and classify of the 1st o'clock phase image, and acquisition is positioned at the locus of the plane closure curve family of cutting apart the rear region border (final curves bunch).

4, the phase high-resolution remote sensing image adopts increment type to cut apart and classify the 2nd to T the time

A) high-resolution remote sensing image of phase still adopts cutting apart and disaggregated model described in the step 2 to the 2nd to T the time, but adopt new starting condition and new image action scope, new starting condition method to set up as shown in Figure 3, its expression: the closed curve family in plane after the phase high-resolution remote sensing image is cut apart during with j is the phase high-resolution remote sensing image preceding initial profile of cutting apart and classify during as the j+1 that is adjacent; For example, finish the 1st o'clock phase high-resolution remote sensing image cut apart and classify (its initialization utilize GIS data auxiliary finish) after, the locus of the closed curve family in plane after will cutting apart is the high-resolution remote sensing image of the phase preceding initial profile of cutting apart and classify during directly as the 2nd, when finishing the 2nd the phase high-resolution remote sensing image cut apart with assorting process after, closed curve family locus, plane after will cutting apart is the high-resolution remote sensing image of the phase preceding initial profile of cutting apart and classify during again as the 3rd, by that analogy, phase image during to after the 3rd o'clock phase high-resolution remote sensing image each, it is cut apart and the preceding initialization mode of classify is all adopted in such a way; Compare with the initial method described in the step 3, this method makes the initial profile locus of each level set function more approach final cutting apart and classification results, and the phase image is cut apart and classified the required time when having reduced j+1.Still adopt the boundary condition described in the step 3, the condition of convergence, discretize mode, and parameter value solution procedure 2 described in cutting apart and disaggregated model, up to reaching new convergence state, the high-resolution remote sensing image of phase cuts apart and classifies when realizing j+1.Still adopt concrete the 2nd step of implementation content c cutting apart with classification results of the high-resolution remote sensing image of phase to the 2nd to T the time) in the expression mode of model solution of elaboration.In above elaboration, T be the change-detection image the time mutually the sum, phase number when its value depends on the image that the change-detection task is used in the practical application, but must satisfy T 〉=2, that is, need at least two not simultaneously the image of phase could use this method to finish the change-detection task, but the time number of phases big more, cumulative errors is big more, causes the change-detection precision to reduce.In addition, because the high-definition remote sensing image data amount is big, and the problem of obtaining difficulty, so this T value can not be too big.Generally get 2～10.In more than setting forth, j is an integer variable, the high-resolution remote sensing image of phase when representing j, and span is for being 1≤j≤T-1.

When b) the classification number of image increases or reduces mutually when two of j and j+1 are adjacent, need increase new level set function on previous level set function classification base plinth or reduce existing certain level set function, the level set function classification number of the level set function initial position of increase and minimizing can be determined by following dual mode: 1) fixed according to priori; 2) by the simple difference of these two images of phase when adjacent, region of variation instructs the setting of the initial profile of newly-increased level set function as the priori of newly-increased classification.

5, the change-detection of phase image adjacent two time

When having finished two of j and j+1 when adjacent mutually after the cutting apart and classify of image, begin the change-detection of this 2 o'clock phase image.For example, during T=3, obtained the 1st when adjacent with the 2 two mutually image cut apart with classification results after, begin these two images are carried out change-detection, the cutting apart and classify of the image of phase when having finished the 3rd, the image that begins the 2nd when adjacent with the 3 two mutually carries out change-detection; Shown in accompanying drawing 4, concrete detection method is: with the classification results of any one image wherein is reference, contrast the classification results of this two phases when adjacent image, if change has taken place in all or part of class label of certain cutting object in the classification results of another image, variation has taken place in the type of ground objects that then shows this object, and the locus and the category attribute of output all changes pixel.

6, returned for the 4th step, the process of circulation execution in step 4 to 6, the phase image cuts apart and classification and change-detection task when having finished all T.

Claims (1)

1. cutting apart and classification and variety detection integration method of a high-resolution remote sensing image, its feature may further comprise the steps:

One, the pre-service of high-resolution remote sensing image:

The phase high-resolution remote sensing image carries out radiant correction respectively during a) to each; And the phase high-resolution remote sensing image carries out registration during to behind the radiant correction each, and then with the GIS data registration of areal;

B), form m dimensional feature vector image to high-resolution remote sensing image texture feature extraction and spectral signature after radiant correction and the registration

Wherein, m is the dimension of the eigenvector image of extraction, and m 〉=1; u _kThe characteristic image of representing k component, 1≤k≤m; Utilize the multiple dimensioned incorgruous diffusion technique of vector total variation to carry out filtering then;

Two, foundation is cut apart and disaggregated model to the eigenvector image after the pre-service:

At first, set up multilevel collection evolution equations, i.e. partial differential equations:

$\left\{\begin{array}{c}\frac{{\&PartialD;\mathrm{\&Phi;}}_1}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_1(x,y)+{\mathrm{\&phi;}}_1(x,y)+{\mathrm{\&mu;Cur}}_{{\mathrm{\&Phi;}}_1}]\left|\right|{\&dtri;\mathrm{\&Phi;}}_1(x,y)\left|\right|\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \frac{{\&PartialD;\mathrm{\&Phi;}}_i}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_i(x,y)+{\mathrm{\&phi;}}_i(x,y)+{\mathrm{\&mu;Cur}}_{{\mathrm{\&Phi;}}_i}]\left|\right|{\&dtri;\mathrm{\&Phi;}}_i(x,y)\left|\right|\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \frac{{\&PartialD;\mathrm{\&Phi;}}_{N-1}}{\&PartialD;t}(x,y)=[-{\mathrm{\&xi;}}_{N-1}(x,y)+{\mathrm{\&phi;}}_{N-1}(x,y)+{\mathrm{\&mu;Cur}}_{{\mathrm{\&Phi;}}_{N-1}}]\left|\right|{\&dtri;\mathrm{\&Phi;}}_{N-1}(x,y)\left|\right|\end{array}\right.---\left(1\right)$

In the formula (1),

Be implicit expression N-1 dimension level set function vector, N is the atural object classification number that high-resolution remote sensing image comprises; Function ξ _i(x y) is i classification C among the Ω of bidimensional image territory _iStatistical nature, computing formula is:

${\mathrm{\&xi;}}_i(x,y)=\frac{1}{m}\underset{k=1}{\overset{m}{\mathrm{\&Sigma;}}}{(u_k(x,y)-c_{\mathrm{ik}})}^4,1\&le;i\&le;N---\left(2\right)$

In the formula (2), c _IkBe i classification C on the k dimensional feature vector image among the Ω of bidimensional image territory _iAverage, computing formula is:

$c_{\mathrm{ik}}=\frac{{\&Integral;}_{C_i}{u}_k(x,y)\mathrm{dxdy}}{{\&Integral;}_{C_i}\mathrm{dxdy}}---\left(3\right)$

In the formula (3), integral sign

The expression limit of integration is i classification C _iThe pixel set that comprises; In the formula (1), function phi _i(x, computing formula y) is:

${\mathrm{\&phi;}}_i(x,y)={\mathrm{\&xi;}}_{i+1}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)>0\}}(x,y)+{\mathrm{\&xi;}}_{i+2}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)\&le;0\}}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+2}(x,y,t)>0\}}(x,y)+\&CenterDot;\&CenterDot;\&CenterDot;$

$+{\mathrm{\&xi;}}_{N-1}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)\&le;0\}}(x,y)\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-2}(x,y,t)\&le;0\}}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-1}(x,y,t)>0\}}(x,y)---\left(4\right)$

$+{\mathrm{\&xi;}}_N(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{i+1}(x,y,t)\&le;0\}}(x,y)\&CenterDot;\&CenterDot;\&CenterDot;{\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-2}(x,y,t)\&le;0\}}(x,y){\mathrm{\&chi;}}_{\{{\mathrm{\&Phi;}}_{N-1}(x,y,t)\&le;0\}}(x,y)$

In the formula (1),

Be the mean curvature of i level set function, computing formula is: ${\mathrm{Cur}}_{{\mathrm{\&Phi;}}_i}=\&dtri;\&CenterDot;\frac{{\&dtri;\mathrm{\&Phi;}}_i}{\left|\right|{\&dtri;\mathrm{\&Phi;}}_i\left|\right|},$

Be gradient operator; Parameter μ is the weights of mean curvature;

Be the indicative function of i level set function, be defined as:

Specifying constraint:

$C_i(x,y,t)=\underset{1\&le;i\&le;N-1}{\min }\left\{{\mathrm{\&xi;}}_i(x,y,t)\right\}---\left(5\right)$

In the formula (5), C _i(x, y, t) in the t time iteration of expression, ξ _i(x, y, t) pixel (x, class label y) of i level set function correspondence of value minimum; When multilevel collection evolutionary process satisfies the condition of convergence, class label C _iVariation remain static;

The partial differential equations (1) that has constraint condition (5) is exactly cutting apart of high-resolution remote sensing image and disaggregated model, and separating in the following ways of this model represented: the 1st is expressed as to the N-1 class: C _i(x, y, t)=(x, y) ∈ Ω | Φ _i(x, y, t)〉0}, and wherein, 1≤i≤N-1, the N class is expressed as: $C_N(x,y,t)=\{(x,y)\&Element;\mathrm{\&Omega;}|\underset{i=1}{\overset{N-1}{\mathrm{}\&cap;\mathrm{}}}{\mathrm{\&Phi;}}_i(x,y,t)\&le;0\},$ And the zone that segmentation result is exactly the zero level collection of all level set functions to surround;

Three, the cutting apart and classify of phase high-resolution remote sensing image the 1st time:

A) utilize method of finite difference to carry out discretize to cutting apart with disaggregated model, boundary condition is the Neumann boundary value condition; Parameter in input type (1)-(3): the atural object classification number N that image comprises gets positive integer, and N 〉=2, and like this, the required level set function number of cutting apart with disaggregated model is N-1; The weights μ of mean curvature=(0.1) p * 255 ², p is an integer, span is 0～5; Spatial sampling interval delta x=1 on the row, column direction, Δ y=1, i.e. a pixel; Iteration time step delta t depends on the evolution speed of spatial sampling interval and level set function, and its value must satisfy CFL condition;

B) in the high-resolution remote sensing image level set of the 1st o'clock phase being cut apart and is classified, cut apart with the starting condition method to set up of disaggregated model as follows: the GIS of this area data behind registration are down auxiliary, the initial profile of each level set function is set, concrete grammar: utilize the locus of each atural object classification in the GIS data and attribute to come initial circle of picture, and with each the initial round starting condition of symbolic distance conversion of carrying out as each level set function, adopt in the step 3 a) described boundary condition and discretize mode and parameter value realize the 1st o'clock mutually the level set of image cut apart and classify, cut apart and obtain cutting apart when classify end after plane closure curve family;

Four, the high-resolution remote sensing image of phase adopts increment type to cut apart and classification the 2nd to T the time: mutually high-resolution remote sensing image still adopts cutting apart and disaggregated model described in the step 2 to the 2nd to T the time, but adopt new starting condition, new starting condition method to set up is as follows: the closed curve family in plane after the phase high-resolution remote sensing image is cut apart during with j is the phase high-resolution remote sensing image preceding initial profile of cutting apart and classify during as the j+1 that is adjacent; Still adopt in the step 3 a) when described boundary condition and discretize mode and parameter value are realized j+1 cutting apart and classifying of mutually high-resolution remote sensing image, the closed curve family in plane after obtaining cutting apart when finishing cut apart and classify; Wherein: T be the change-detection image the time mutually the sum; J is an integer variable, the high-resolution remote sensing image of phase when representing j, and span is for being 1≤j≤T-1;

Five, the change-detection of high-resolution remote sensing image mutually when two of j and j+1 are adjacent: from j=1, when j during with j+1 mutually high-resolution remote sensing image cut apart finish with assorting process after, begin to these two that the phase high-resolution remote sensing image carries out change-detection when adjacent, concrete grammar is: the classification results that contrasts these two high-resolution remote sensing images of phase when adjacent, with wherein any one the time phase high-resolution remote sensing image classification results be reference, another the time phase high-resolution remote sensing image classification results in, with the object after cutting apart is comparing unit, if change has taken place in the whole pixels of any one object or the class label of part pixel, variation has taken place in the type of ground objects that then shows this object, and the locus and the category attribute of output all changes pixel;

Six, returned for the 4th step, the process of circulation execution in step four to six, the phase image cuts apart and classification and change-detection task when having finished all T.

Images