CN100351647C - Area feature variation detection method based on remote sensing image and GIS data - Google Patents

Abstract

The present invention relates to an area ground object variation detection method based on a remote sensing image and GIS data, which comprises the steps: step 1, pretreating the GIS data; step 2, pretreating the remote sensing image; step 3, calculating the label point of an area ground object; step 4, predicting the matching point of the remote sensing image, which corresponds to the label point of the area ground object; step 5, defining a group of similar characteristic amounts from thickness to thinness; step 6, extracting the matching area ground object on the image according to the iteration of the similar characteristic amounts; step 7, judging whether change happens in the ground object or not according to the threshold value of the predefined similar characteristic amount, marking the ground object if the change happens, or calculating the polygonal barycenter coordinate of the ground object as control points, and returning the step 4. The method is used to operate the GIS data of the interesting ground object to detect the variation of the ground object and obtain the control points to the ground object without variation, and the control points can be used to realize the automatic registration of the image and the GIS data.

Description

Area feature variation detection method based on remote sensing image and GIS data

Technical field

The invention belongs to the remote sensing image application that combines with GIS (Geographic Information System), relate to the whole iterative solution of a kind of area feature change-detection and ask method based on remote sensing image and GIS data.

Background technology

Natural variation and human comings and goings are all changing face of land view and are utilizing form every day.The quick growth and the development of urbanization of population have been quickened this variation.These variations will produce far-reaching influence to earth resource and environment; therefore monitoring the face of land in time, effectively changes; upgrading relevant Geographic Information System, is very necessary for functional departments such as resource management and planning and environmental protection provide the foundation of science decision.For understanding the satellite earth observation technology that the earth grows up, monitor the best-of-breed technology means that the face of land changes beyond doubt.Modern Remote Sensing Technical has entered a new stage that multiple earth observation mass data can fast, in time be provided.Will be from the remote sensing image of number in TB of present reception every day, sort out our interested data, rely on traditional artificial decomposition method obviously not all right, this just need guide our notice into our interested place by allowing computing machine understand image someway and detecting where variation (change-detection) has taken place.After the spatial data framework builds up, in order to realize the data fast updating, keep its trend of the times, it is particularly important that change-detection seems, and change-detection is one of remote sensing and GIS field emphasis of future research (Li Deren 2003) automatically.

Existing change detecting method mainly still rests on the data pilot of Pixel-level, lacks the feature level change detecting method (Li Deren 2003) of knowledge elicitation.At present the change detecting method model mainly be to not simultaneously the remote sensing image of phase handle, its precondition is the significant change that the time phase change of atural object can cause pixel value on the image.This class change detecting method has only utilized gradation of image information, does not utilize geological information and other knowledge of detected object.Because gradation of image has certain uncertainty to the expression of atural object, thereby make the result of change-detection a lot of pseudo-variations occur.Although there have been many change detecting methods, select a suitable method still very difficult for a certain specific application and survey region.Comprised abundant semanteme and non-semantic information in the GIS data, integrated remote sensing images and GIS data are carried out one of the trend (D.Lu 2004) that the change-detection analysis is change detecting method development in recent years.

General change detecting method step can reduce at present: data pre-service, data registration, feature extraction (based on the change detecting method of feature level), change-detection, interpretation of result, these steps are independently to carry out, the error in each step all can be delivered to next step, wherein data registration and feature extraction have critical influence to change-detection, make that the precision as a result and the reliability of change-detection are uncontrollable, and testing process also relatively blindly.The present invention proposes a kind of change detecting method based on the feature level for this reason, the committed step that influences the change-detection result: the synchronous iterative of data registration, feature extraction and change-detection, constraint mutually, revise.Utilize knowledge channeling conduct in the GIS data to reduce the blindness of feature extraction and change-detection, improve result's reliability and precision.

Summary of the invention

The purpose of this invention is to provide a kind of area feature variation detection method based on remote sensing image and GIS data, this method can realize multi-source, multisensor, not simultaneously the number of phases according between change-detection.

Technical scheme provided by the invention is that a kind of area feature variation detection method based on remote sensing image and GIS data may further comprise the steps:

One, to the pre-service of GIS data;

I. from the GIS data, choose area feature, and extract the polygon of expressing its geometric configuration by actual one-tenth figure ratio expression;

II. calculate each polygonal following similar features amount among the I: the area of the minimum boundary rectangle of polygon; The minimum boundary rectangle of polygon wide and high; Polygonal area; Polygonal girth; Polygonal shape coding;

III. determine polygonal Label point (point promptly, down with) among the I by following step:

1) carries out rasterizing according to each polygonal coordinate, generate the image of polygon boundary rectangle size;

2) with erosion operation (Erosion) to 1) in image carry out corrosion treatment to convex polygon image and become straight line; The concave polygon figure becomes the curve of one or more UNICOM;

3) to 2) result that generates, from these straight lines or curve, choose an intermediate point, as the Label point; Perhaps continue to utilize erosion operation, the line image that obtains is carried out computing, become a point, at last this that obtains a Label point as image up to line;

Two, to the remote sensing image pre-service;

Treatment step is as follows:

1) initialization i=0, input original image f _i(x, y);

2) with two-dimensional convolution operator and image f _i(x y) carries out convolution, obtains f _I+1(x, y):

3) calculate first wavelet coefficient: w _I+1(x, y)=f _i(x, y)-f _I+1(x, y);

4) if i＜n, wherein n is given decomposition number of times, i=i+1 returns 2);

5) repeat two .2), two .3), two .4) until i=n;

6) choose 1-3 small echo face as remote sensing image pre-service result;

The method of mirror image symmetry is adopted in the processing on remote sensing image border, that is:

Row to: f (i, j)=f (i, j);

F (i+k, j)=f (i-k, j) i＜=N wherein, k=1,2 ..., N is total line number of image; Row to: f (i ,-j)=f (i, j);

F (i, j+k)=f (i, j-k) j＜=N wherein, k=1,2 ..., N is total columns of image;

Three, set up rough geometry site:,, set up the rough geometry site between GIS data and the remote sensing image by the geometric transformation model parameter between least square method calculating remote sensing image and the GIS data according to the reference mark of determining;

Four, read through the pretreated polygon of GIS data, according to the geometry site of setting up in three, the position of prediction GIS polygon Label point corresponding same place on remote sensing image, with this point is seed points, takes the area feature polygon feature on the adaptive iteration region growing method extraction remote sensing image; The threshold value of gray consistency is started from scratch along with the iterations increase progressively increases (increasing the extracted region threshold value), polygon to each extraction calculates minimum boundary rectangle area (analog quantity), difference between the minimum boundary rectangle area of polygon that extracts on the calculating remote sensing image and GIS polygon is up to satisfying preselected threshold condition; Note the polygon that each satisfies threshold value, as polygonal candidate polygon of the same name among the GIS;

Five, the change-detection of area feature:

1. calculated candidate polygonal similar features amount of the same name: the minimum boundary rectangle of polygon wide and high; Polygonal area; Polygonal girth; Polygonal shape coding;

2. wide and height is the similar features amount with the minimum boundary rectangle of polygon, selects the polygon of feature similarity from five, 1 result;

3. be similarity measure with polygonal area, from five, 2 results, select the polygon of feature similarity;

4. be similarity measure with polygonal girth, from five, 3 results, select the polygon of feature similarity;

5. be similarity measure with polygonal shape coding, from five, 4 results, select the polygon of feature similarity; The polygon of the similarity maximum of selected shape coding is a feature of the same name;

6. if the result in five, 5 does not satisfy pre-set threshold, just think that variation has taken place this area feature, this GIS polygon of mark is a region of variation; Otherwise (maximum similar value is greater than the similarity threshold value) this area feature does not change, and the of the same name polygonal center of gravity of getting on this polygon and the remote sensing image is the reference mark, participates in next area feature change-detection, optimizes rough geometry site;

Six, each the GIS polygon to be detected that extracts among its rapid I is carried out iterative processing according to step 3 to five, all mate until all polygons and finish.At last, the result who detects is carried out precision evaluation.

The present invention carries out aforesaid operations to the GIS data of atural object interested, can detect variation atural object, and the atural object for not having to change can obtain the reference mark, utilizes these reference mark can realize the autoregistration of image and GIS data.The present invention utilizes the knowledge channeling conduct in the GIS data to reduce the blindness of feature extraction and change-detection, improves result's reliability and precision.The present invention can realize multi-source, multisensor, not simultaneously the number of phases according between change-detection.

Description of drawings

Accompanying drawing is a schematic flow sheet of the present invention.

Embodiment

Referring to accompanying drawing, the present invention includes following steps:

One, to the pre-service of GIS data

I. from the GIS data, choose area feature, and extract the polygon of expressing its geometric configuration by actual one-tenth figure ratio expression.

II. sort by area is descending according to the polygon coordinate reference area, and to all polygons.

III. calculate the polygonal similar features amount that each group is used to mate.Step 4 is stated in the calculating of similarity characteristic quantity as follows.

IV. determine polygonal " Label " point.Each polygon according to following steps, is calculated the Label point:

1) in calculator memory, carries out rasterizing according to each polygonal coordinate, generate the bianry image of polygon boundary rectangle size.

2) with the erosion operation in the mathematical morphology (Erosion) to 1) in the result carry out repeatedly corrosion treatment.In the process of corroding, generally getting the corrosion structure element is [1 1; 1 1].The symmetry of corrosion structure element mainly is the direction of control corrosion.

3) if image is a convex polygon, in process step 2) after, can become straight line: if figure is a concave polygon, the result of corrosion can be the curve of one or more UNICOM.

4) to 3) result that generates, when choosing the Label point, two kinds of methods are arranged: the one, from these straight lines or curve, choose an intermediate point, as the Label point: the 2nd, continue to utilize erosion operation, the line image that obtains is carried out computing, become a point up to line, this point that obtains at last just can be used as the Label point of image.To 3) in result when corroding, the corrosion structure element is commonly defined as [1 0; 0 1].

Two, based on

The remote sensing image pre-service of wavelet decomposition

Ask the polygon feature that requires to extract area feature in the method based on the global solution of remote sensing image and GIS change-detection, wish that for this reason the image preprocessing process strengthens marginal information, level and smooth homogeneous region when restraining picture noise from image.Here it is a kind of to have used Research on wavelet transform method The remote sensing image preprocess method of wavelet decomposition.The method principle is as follows:

If function ψ (x) ∈ is L ²(R) (L ²(R) be the square integrable space), and ψ (x) satisfies

$C_{\mathrm{\&psi;}}=2\mathrm{\&pi;}\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}\frac{{\left|\stackrel{\&OverBar;}{\mathrm{\&psi;}}\left(\mathrm{\&omega;}\right)\right|}^2}{\left|\mathrm{\&omega;}\right|}\mathrm{d\&omega;}<+\&infin;---(2-1)$

Or $\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}\mathrm{\&psi;}\left(x\right)\mathrm{dx}=0---(2-2)$

In the formula: ψ (ω) is the Fourier transform of ψ (x).

When ψ (x) satisfies (2-1) or (2-2) formula, and can faster restrain the time, we claim that ψ (x) is basic small echo.After ψ (x) is through flexible a and translation b operation:

${\mathrm{\&psi;}}_{a,b}\left(x\right)={\left|a\right|}^{-\frac{1}{2}}\mathrm{\&psi;}\left(\frac{x-b}{a}\right)---(2-3)$

In the formula: a, b ∈ R and a ≠ 0

At this moment, ψ a, b (x) is called small echo.

The wavelet transformation of function f (x) is defined as:

$\mathrm{wf}(a,b)=\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}f\left(x\right)\&CenterDot;{\left|a\right|}^{\frac{1}{2}}\stackrel{\&OverBar;}{\mathrm{\&psi;}}\left[a^{-1}(x-b)\right]\mathrm{dx}---(2-4)$

In the formula: f (x) ∈ L ²(R), b ∈ R, ψ (x) is ψ, complex conjugate (x).

Introduced the definition of continuous wavelet transform above, continuous wavelet can not be used for programming and calculate, must discretize.The method of the discretize of continuous wavelet transform has multiple, and is still, all effective for the not every discrete method of different problems.Famous Mallat algorithm is to utilize orthogonal basis to carry out tower decomposition, but through the size of the image after this method conversion variation has taken place, and this friendshipization is disadvantageous often in some image processing process, for example: pattern-recognition, multi-source visual fusion etc.

In order to make image constant, adopt a kind of being called as here through size behind the wavelet transformation

The method of wavelet decomposition becomes different small echo planes to picture breakdown.

The basic thought of wavelet algorithm is to be signal or picture breakdown approximate signal on the different frequency passage and the detail signal under each yardstick.This detail signal is called little corrugated, and its image size is identical with original image size.

For one-dimensional signal C (x), suppose { C ₀(x) } be the scalar product of signal C (x) and scaling function φ (x), scaling function is actually a low-pass filter.Obtain C after the filtering for the first time of signal C (x) process ₁(x), w ₁(x)=C ₀(x)-C ₁(x) comprise information between these two yardsticks, w ₁(x) being called the first little corrugated, also is the result of the wavelet transformation of corresponding scaling function.And wavelet function ψ (x) has following relation with scaling function φ (x):

$\frac{1}{2}\mathrm{\&psi;}\left(\frac{x}{2}\right)=\mathrm{\&phi;}\left(x\right)-\frac{1}{2}\mathrm{\&phi;}\left(\frac{x}{2}\right)---(2-5)$

Differ twice between the adjacent yardstick, obtain C through after i the filtering _i(x) be:

$C_i\left(x\right)=\underset{1}{\mathrm{\&Sigma;}}h\left(1\right){\&CenterDot;C}_{i-1}(x+2^{i-1}1)---(2-6)$

The wavelet transform wavelet coefficient is:

w _i(x)＝C _i-1(x)-C _i(x) (2-7)

w _i(x) be the wavelet coefficient under the yardstick i (little corrugated), C _i(x) be the approximate signal under the i yardstick, h is a low-pass filter, establishes an equation under it and scaling function φ (x) satisfy:

$\frac{1}{2}\mathrm{\&phi;}\left(\frac{x}{2}\right)=\underset{1}{\mathrm{\&Sigma;}}h\left(1\right)\&CenterDot;\mathrm{\&phi;}(x-1)---(2-8)$

Discrete wavelet decomposes signal, generates little corrugated { wi} and the approximate signal sum of one group of adjacent different resolution.

If the scaling function of disclosing in the selection linearity is promptly:

φ (x) if=1-|x| x ∈ [1,1]

φ (x) if=0 x  [1,1] (2-9)

So can calculate h (1)=1/4, h (0)=1/2, h (1)=1/4:

$\frac{1}{2}\mathrm{\&phi;}\left(\frac{x}{2}\right)=\frac{1}{4}\mathrm{\&phi;}(x+1)+\frac{1}{2}\mathrm{\&phi;}\left(x\right)+\frac{1}{4}\mathrm{\&phi;}(x-1)$

$C_{i+1}\left(x\right)=\frac{1}{4}{C}_i(x-2^i)+\frac{1}{2}{C}_i\left(x\right)+\frac{1}{4}{C}_i(x+2^i)---(2-10)$

Above-mentioned

Thereby discrete wavelet decomposes and to be easy to be generalized to 3 * 3 the convolution operator that two dimension obtains a two dimension:

$\left[\begin{array}{ccc}\frac{1}{16}& \frac{1}{8}& \frac{1}{16}\\ \frac{1}{8}& \frac{1}{4}& \frac{1}{8}\\ \frac{1}{16}& \frac{1}{8}& \frac{1}{16}\end{array}\right]$

If adopt the spline scale function B3 time, the two-dimensional convolution operator is so:

$\left[\begin{array}{ccccc}\frac{1}{256}& \frac{1}{64}& \frac{3}{128}& \frac{1}{64}& \frac{1}{256}\\ \frac{1}{64}& \frac{1}{16}& \frac{3}{32}& \frac{1}{16}& \frac{1}{64}\\ \frac{3}{128}& \frac{3}{32}& \frac{9}{64}& \frac{3}{32}& \frac{3}{128}\\ \frac{1}{64}& \frac{1}{16}& \frac{3}{32}& \frac{1}{16}& \frac{1}{64}\\ \frac{1}{256}& \frac{1}{64}& \frac{3}{128}& \frac{1}{64}& \frac{1}{256}\end{array}\right]$

Performing step is as follows:

1) initialization i=0, input original image f _i(x, y);

2) (x is y) with image f with wave filter h _i(x y) carries out convolution, obtains f _I+1(x, y):

f _i+1(x，y)＝f _i(x，y)×h(x，y)：

3) carry out the small echo variation first time, obtain first wavelet coefficient:

w _i+1(x，y)＝f _i(x，y)-f _i+1(x，y)；

4) if i＜n (n is given decomposition number of times), i=i+1 returns 2);

5) repeat 2), 3), 4) and until i=n.

The method of mirror image symmetry is adopted in the processing on border, that is:

Row to: f (i, j)=f (i, j):

F (i+k, j)=f (i-k, j) i＜=N wherein, k=1,2 ..., N is total line number of image;

Row to: f (i ,-j)=f (i, j):

F (i, j+k)=f (i, j-k) j＜=N wherein, k=1,2 ..., N is total columns of image.

Select suitable small echo face number in actual applications as required for use.

Three, geometric transformation CALCULATION OF PARAMETERS

Root pick reference mark is by the geometric transformation model parameter between least square method calculating remote sensing image and the GIS data, by the geometry site between geometric transformation modelling GIS data coordinates and the remote sensing image pixel.Geometric transformation model commonly used has: affined transformation, polynomial transformation and perspective projection transformation (collinearity equation) etc.

The geometric transformation model of image can adopt three kinds of modes respectively according to different actual conditions: (1) is for the more smooth area of landform, adopt earlier simple polynomial transformation as the rough model of remotely sensing image geometric conversion, treat to have determined that by mating using strict image rectification model again behind abundant, the high-precision reference mark carries out high precision differential rectify; (2) adopt the strictness of remotely sensing image geometric distortion to correct model at the very start, along with the raising of reference mark number and the precision model parameter of constantly refining; (3) mixed type of employing (1) and (2), adopt simple polynomial transformation as the rough model of remotely sensing image geometric conversion earlier, treat by coupling got access to can separate ask the minimum reference mark of strict model parameter number after, adopt the strict correction model of image distortion, again along with the increase of reference mark number is carried out iterative refinement to the image rectification model parameter of strictness.

Four, similar features amount

The definition of similar features amount is very crucial, and it is a yardstick of judging feature of the same name.The definition of similar features amount is directly connected to reliability, stability and the uniqueness of testing result, has also determined the size of calculated amount.

Here define following five similar features amounts for the polygon feature of area feature and describe its similarity:

1, the area of the minimum boundary rectangle of polygon;

2, the minimum boundary rectangle of polygon is wide and high;

3, polygonal area;

4, polygonal girth;

5, polygonal shape coding;

These feature similarity amounts have been formed a similarity measure set from coarse to fine together, though list can not uniquely be determined feature of the same name with regard to each amount of estimating, but the similar mensuration of this group can narrow down to candidate feature in the very little scope, add the constraint of the rough geometric position of calculating in three, just can uniquely determine feature of the same name.

Five, area feature polygon Feature Extraction

Read through the pretreated polygon of GIS data, according to the geometry site of setting up in three, the position of prediction GIS polygon Label point corresponding same place on remote sensing image, with this point is seed points, takes the area feature polygon feature on the adaptive iteration region growing method extraction remote sensing image; The threshold value of gray consistency is started from scratch along with the iterations increase progressively increases (increasing the extracted region threshold value), polygon to each extraction calculates minimum boundary rectangle area (analog quantity), difference between the minimum boundary rectangle area of polygon that extracts on the calculating remote sensing image and GIS polygon is up to satisfying preselected threshold condition; Note the polygon that each satisfies threshold value, as polygonal candidate polygon of the same name among the GIS.So just realized the image-region polygon feature extracted in self-adaptive under the GIS data pilot automatically.

Six, search strategy

The geometric transformation result of ordering with GIS polygon Label is guiding, in order to take into account counting yield and accuracy, adopts following search strategy:

1, progressively mates from big to small by area of a polygon.

The geometric error that the polygon that area is big can allow is big, by the descending precision that progressively improves the geometric transformation model of area, can guarantee to drop in the image-region of the same name through the GIS polygon Label point that the image geometry transformation model calculates.

2, the coupling to feature of the same name adopts layering and matching, the progressively strategy of refinement.Concrete grammar is as follows:

1) the geometric transformation result of ordering with the polygonal Label of GIS is the seed points that image polygonal region of the same name extracts, and is similarity measure with the minimum boundary rectangle area of GIS polygon, progressively increases the gray consistency threshold value, extracts the image polygonal region.The similar tolerance threshold value of minimum boundary rectangle area can preestablish, and the setting of threshold value can be looser.Extract, note the polygon that each satisfies threshold value, from image as candidate's polygon feature of the same name.

2) wide and height is a similarity measure with the minimum boundary rectangle of polygon, from 1) select the polygon of feature similarity the result.

3) be similarity measure with polygonal area, from 2) select the polygon of feature similarity the result.

4) be similarity measure with polygonal girth, from 3) select the polygon of feature similarity the result.

5) be similarity measure with polygonal shape coding, from 4) select the polygon of feature similarity the result.The polygon of the similarity maximum of selected shape coding is a feature of the same name.

Seven, the change-detection of area feature

If the result 6 5) does not satisfy predefined similarity threshold value, just think that variation has taken place this area feature, this GIS polygon of mark; Otherwise (being that maximum similar value is greater than the similarity threshold value) this area feature does not change, the of the same name polygonal center of gravity of getting on this polygon and the remote sensing image is the reference mark, participate in next area feature change-detection, optimize rough geometry site, the GIS polygon to be detected that extracts in step 1 I all disposes.At last, the result who detects is carried out precision evaluation.

When the present invention is used for not doing change-detection between the GIS vector data of phase simultaneously, can omits for second, five steps, and directly utilize phasor coordinate to describe polygon atural object.

When the present invention was used for not doing change-detection between the raster data of phase simultaneously, the raster data of phase carried out atural object and extracts during earlier to one of them, handles by method of the present invention again.

Claims (3)

1. the area feature variation detection method based on remote sensing image and GIS data is characterized in that: may further comprise the steps

One, to the pre-service of GIS data;

I. from the GIS data, choose area feature, and extract the polygon of expressing its geometric configuration by actual one-tenth figure ratio expression;

II. calculate each polygonal following similar features amount among the I: the area of the minimum boundary rectangle of polygon; The minimum boundary rectangle of polygon wide and high; Polygonal area; Polygonal girth; Polygonal shape coding;

III. determine point in polygonal among the I by following step:

1) carries out rasterizing according to each polygonal coordinate, generate the image of polygon boundary rectangle size;

2) with erosion operation to 1) in image carry out corrosion treatment to convex polygon image and become straight line; The concave polygon figure becomes the curve of one or more UNICOM;

3) to 2) result that generates, from these straight lines or curve, choose an intermediate point, as interior point; Perhaps continue to utilize erosion operation, the line image that obtains is carried out computing, become a point, at last this that obtains an interior point as image up to line;

Two, to the remote sensing image pre-service,

Treatment step is as follows:

1) initialization i=0, input original image f _i(x, y);

2) with two-dimensional convolution operator and image f _i(x y) carries out convolution, obtains f _I+1(x, y);

Wherein the two-dimensional convolution operator is: $\left[\begin{array}{ccccc}\frac{1}{256}& \frac{1}{64}& \frac{3}{128}& \frac{1}{64}& \frac{1}{256}\\ \frac{1}{64}& \frac{1}{16}& \frac{3}{32}& \frac{1}{16}& \frac{1}{64}\\ \frac{3}{128}& \frac{3}{32}& \frac{9}{64}& \frac{3}{32}& \frac{3}{128}\\ \frac{1}{64}& \frac{1}{16}& \frac{3}{32}& \frac{1}{16}& \frac{1}{64}\\ \frac{1}{256}& \frac{1}{64}& \frac{3}{128}& \frac{1}{64}& \frac{1}{256}\end{array}\right]$

3) calculate first wavelet coefficient: w _I+1(x, y)=f _i(x, y)-f _I+1(x, y);

4) if i＜n, wherein n is given decomposition number of times, i=i+1 returns two .2);

5) repeat two .2), two .3), two .4) until i=n;

6) choose 1-3 small echo face as remote sensing image pre-service result;

The method of mirror image symmetry is adopted in the processing on remote sensing image border, that is:

Row to: f (i, j)=f (i, j);

F (i+k, j)=f (i-k, j) i＜=N wherein, k=1,2 ..., N is total line number of image;

Row to: f (i ,-j)=f (i, j);

F (i, j+k)=f (i, j-k) j＜=N wherein, k=1,2 ..., N is total columns of image;

Three, according to the reference mark of determining, by the geometric transformation model parameter between least square method calculating remote sensing image and the GIS data;

Four, take the adaptive iteration region growing method to extract area feature polygon feature on the remote sensing image; The seed points of region growing is provided by the corresponding geometric transformation result of the polygonal interior point of GIS on the remote sensing image, the threshold value of gray consistency is started from scratch along with the iterations increase progressively increases, polygon to each extraction calculates minimum boundary rectangle area, difference between the minimum boundary rectangle area of polygon that extracts on the calculating remote sensing image and GIS polygon is up to satisfying preselected threshold condition; Note the polygon that each satisfies threshold value, as polygonal candidate polygon of the same name among the GIS;

Five, the change-detection of area feature:

1) calculated candidate polygonal similar features amount of the same name: the minimum boundary rectangle of polygon wide and high; Polygonal area; Polygonal girth; Polygonal shape coding;

2) wide and height is the similar features amount with the minimum boundary rectangle of polygon, from five .1) select the polygon of feature similarity the result;

3) be the similarity characteristic quantity with polygonal area, from five .2) select the polygon of feature similarity the result;

4) be the similarity characteristic quantity with polygonal girth, from five .3) select the polygon of feature similarity the result;

5) be the similarity characteristic quantity with polygonal shape coding, from five .4) select the polygon of feature similarity the result; The polygon of the similarity maximum of selected shape coding is a feature of the same name;

6) if five .5) in the result do not satisfy pre-set threshold, just think that variation has taken place this area feature; Otherwise this area feature does not change, and the of the same name polygonal center of gravity of getting on this polygon and the remote sensing image is the reference mark, participates in next area feature change-detection;

Six, each the GIS polygon to be detected that extracts among the step 1 I is carried out iterative processing according to step 3 to five.

2. method according to claim 1 is characterized in that: at step 1 .III.2) in the process of corroding, getting the corrosion structure element is [1 1; 1 1].

3. method according to claim 1 and 2 is characterized in that: at step 1 .III.3) in the process of corroding, getting the corrosion structure element is [1 0; 0 1].

Images