CN103218783A - Fast geometric correction method for satellite remote sensing image and based on control point image database - Google Patents

Abstract

A fast geometric correction method for a satellite remote sensing image and based on a control point image database comprises the following steps: the geographical scope of an image to be corrected is determined; all control points which meet requirements are retrieved in the control point image database according to the determined geographical scope of the image to be corrected; automatic matching and point selection are conducted, and control point pairs used for geometric correction are determined; whether the number and the distribution of the matched control point pairs meet the requirements of the geometric correction is judged, if the judging result is yes, the next step is conducted, and if the judging result is no, the retrieval step is conducted again; a geometric correction model is selected according to the type and imaging features of a sensor; a transformational relation between a pixel coordinate and a geodetic coordinate is established based on the matched control point pairs and the selected geometric correction model; and the geometric correction is conducted based on a digital elevation model and the transformational relation so as to obtain a corrected digital orthoimage. According to the fast geometric correction method for the satellite remote sensing image and based on the control point image database, the automatic image matching technique is introduced into automatic selection of the control points, and automatic geometric correction processing of the satellite remote sensing image can be realized.

Description

Satellite remote sensing images fast geometric correcting method based on the reference mark image database

Technical field

The present invention relates to a kind of remote sensing images correcting method, more specifically relate to a kind of satellite remote sensing images fast geometric correcting method based on the reference mark image database.

Background technology

Along with the development of remote sensing technology, the particularly continuous development of remote sensor technology, more and more wider by the purposes of remote sensing image that remote sensing technology obtained or data.At present, the range of application of remotely-sensed data has expanded to social information's service field, for example, be widely used in aspects such as mapping, agricultural, forestry, geological and mineral, the hydrology and water resource, environmental monitoring, disaster, regional analysis and planning, military affairs, soil utilization.Remote sensing image with accurate geographic coding can provide characters of ground object and the information that needs separately for different fields such as soil, planning, environmental protection, agricultural, forestry, oceans.

Via satellite or flying platform such as aviation platform when obtaining remote sensing image data or other data, the influence of external factor such as can be subjected to weather, daylight, block, simultaneously, the height of flying platform, attitude can change when data acquisition, therefore, carrying out tending to cause problems such as image translation, rotation, convergent-divergent when remote sensing images are taken.In addition, according to optical imaging concept, according to the imaging of central projection mode, therefore ground height rises and falls and will cause the existence of height displacement when imaging during camera imaging.Above-mentioned combined factors can cause the error of remote sensing image, for example droop error, projection error etc.Therefore, before using these remote sensing image/data, need the original remote sensing image that is obtained is carried out orthorectify.

Traditional remote sensing image orthorectify generally comprises: at first by the field operation measurement or from existing topomap data collection ground control point (Ground Control Point, GCP) and digital elevation model (Digital Elevation Model, DEM) information; Then these information are imported the remote sensing or the digital Photogrammetric System of specialty; Next by the system that has loaded GCP and DEM information remote sensing image is carried out orthorectify.

Ground control point GCP (Ground Control Points) is that satellite remote sensing image geometry is corrected and important reference data sources during geo-location.Correct in the processing procedure at remotely sensing image geometric, for reaching certain correction precision, the ground control point of some is absolutely necessary, object coordinates by the reference mark and corresponding picpointed coordinate make up video imaging model and solving model parameter or existing imaging model is optimized compensation finds the solution compensating parameter, finally sets up the correct transformational relation of object space and picture side in the imaging process.

The general pattern that adopts full field survey in reference mark in the traditional work process, need through a series of complicated step such as " collect existing control data, survey that the district is made an on-the-spot survey, reconnaissance is laying of markstone, field operation is measured, interior industry put in order ", though development along with advanced measuring techniques such as GPS-RTK, workload and complexity that field operation is measured significantly reduce, but necessary control point field survey work is still inevitable.

In addition, (with No. three satellites of resource is example because the satellite remote-sensing image coverage is big, coverage is 50 kilometers * 50 kilometers), seek out equally distributed ground control point, usually need testing in the scope of hundreds of kilometer and even thousands of square kilometres, its field process amount is big, labor intensive and material resources are self-evident.In addition, in the multiple zone of disasteies such as earthquake, flood, rubble flow or uninhabited areas such as virgin forests, marsh, desert, survey crew usually can't enter on the spot and measure.

Even the reference mark is by manually choosing remote sensing image and the topomap same place obtains, also existing inefficiency, precision to be difficult to problems such as guarantee.And to areal not simultaneously mutually or the different sensors image correct, can duplicate the situation of reconnaissance again.

Summary of the invention

According to embodiments of the invention, provide a kind of satellite remote sensing images fast geometric correcting method based on the reference mark image database.This method comprises: the geographic range of determining to wait to correct image; At determined geographic range of waiting to correct image, all satisfactory reference mark of retrieval in the image database of reference mark; Automatically mate reconnaissance, the reference mark that is identified for geometric correction is right; Judge the quantity that the reference mark of described coupling is right and distribute whether meeting geometric is corrected requirement, if, then enter next procedure, if not, described searching step then returned; According to sensor type and imaging characteristics, select corresponding geometric correction model; Based on the reference mark of described coupling to the geometric correction modelling pixel coordinate of described selection and the transformational relation of terrestrial coordinate; Carry out geometric correction based on the transformational relation between digital elevation model and pixel coordinate and the terrestrial coordinate, obtain digital orthoimage through orthorectify.

According to the remote sensing images correcting method of the embodiment of the invention, the automatic matching technique of image is incorporated in the choosing automatically of reference mark, can realize that the automatic geometric of remote sensing image is corrected processing.

Description of drawings

In order to be illustrated more clearly in the technical scheme of the embodiment of the invention, will do to introduce simply to the accompanying drawing of embodiment below, apparently, the accompanying drawing in describing below only relates to some embodiments of the present invention, but not limitation of the present invention.

Fig. 1 is the process flow diagram according to the satellite remote sensing images fast geometric correcting method based on the reference mark image database of the present invention;

Fig. 2 a-2c is the exemplary reference mark image sheet and the comparison diagram of No. three satellite images of resource;

Fig. 3 is a kind of indicative flowchart of reference mark image sheet retrieval;

Fig. 4 is a kind of indicative flowchart of reference mark automatic matching method;

Fig. 5 shows one group of example of the bidimensional image under the different scale space;

Fig. 6 shows one group of example of difference of Gaussian metric space (DoG) image;

Fig. 7 is the synoptic diagram that DoG metric space local extremum detects;

Fig. 8 is the synoptic diagram by key point neighborhood gradient information generating feature vector;

Fig. 9 shows the one group example of different images by key point neighborhood gradient information generating feature vector plot;

Figure 10 is based on the synoptic diagram of flow process of the matching algorithm of SIFT descriptor;

Figure 11 is the synoptic diagram that has comprised that the mistake match point is rejected the matching algorithm flow process of step;

Figure 12 shows the imaging mode of linear array push-broom type imaging sensor;

Figure 13 shows picture point and the topocentric mathematical relation in the geometric correction;

Figure 14 shows the coordinate system that strict imaging model may be used;

Figure 15 is the synoptic diagram of the intersection of imaging direction and earth ellipsoid.

Embodiment

For the purpose, technical scheme and the advantage that make the embodiment of the invention is clearer,, the technical scheme of the embodiment of the invention is clearly and completely described below in conjunction with the accompanying drawing of the embodiment of the invention.Obviously, described embodiment is a part of embodiment of the present invention, rather than whole embodiment.Based on described embodiments of the invention, those of ordinary skills belong to the scope of protection of the invention at the every other embodiment that need not to be obtained under the prerequisite of creative work.

Unless define in addition, technical term of Shi Yonging or scientific terminology should be the ordinary meaning that the personage understood that has general technical ability in the affiliated field of the present invention herein.Any order, quantity or importance do not represented in " first " of using in patent application specification of the present invention and claims, " second " and similar word, and just are used for distinguishing different ingredients.Equally, restricted number do not represented in " one " or similar words such as " one " yet, but there is at least one in expression.

Satellite remote sensing images fast geometric correcting method according to the embodiment of the invention based on the reference mark image database, be the basis that is established as, by retrieval, the reference mark of reference mark image database mating automatically and target image is corrected and realized with the reference mark image database.Described reference mark image database is the set of reference mark image sheet.And reference mark image sheet is the expansion at traditional reference mark, that is to say, replace single reference mark as the data of database unit with orthography DOM piece with geographic coordinate information and corresponding digital elevation model DEM piece, thereby generate the reference mark image database.

Fig. 1 schematically shows the process flow diagram according to the remote sensing images geometric correction method of the embodiment of the invention.At step S1, determine to wait to correct the geographic range of image; At step S3,, in the image database of reference mark, retrieve all satisfactory reference mark at determined geographic range of waiting to correct image; At step S5, mate reconnaissance automatically, the reference mark that is identified for geometric correction is right; At step S7, judge the quantity that the reference mark obtained is right and distribute whether meeting geometric is corrected requirement, if, then enter next procedure, if not, then return step S3 and proceed retrieval; At step S9,, select corresponding geometric correction model classification according to sensor type and imaging characteristics; Step S11 is based on the transformational relation of reference mark to, geometric correction modelling pixel coordinate and terrestrial coordinate; At step S13, the transformational relation based between digital elevation model (DEM) and pixel coordinate and the terrestrial coordinate carries out geometric correction to raw video, obtains the digital orthoimage through orthorectify.

Alternatively, before carrying out step S1, image original to be corrected is carried out image slightly correct.This thick correction can be adopted non-conventional method based on the reference mark image database, for example, adopt four the angle point terrestrial coordinates and the image coordinate that write down in the metadata information of raw video, make up one time multinomial model, finish the geometric correction of summary, obtain the correction image that does not correct through the reference mark.By the effective scope of the raw video of PREDICTIVE CONTROL point image sheet correspondence of thick correction, thereby dwindle the hunting zone of image coupling, reduce match time.Therefore, adopt the thick correction of this image can improve the speed of geometric correction of the present invention.

According to the above-mentioned matching algorithm of the embodiment of the invention, in certain error range, can realize rapid automatized reconnaissance operation.Alternatively, in step S5, adopt artificial assisting sifting to combine, can improve the selection quality at reference mark, and then improve the precision of geometric correction with automatic coupling.

In addition, alternatively, after step S13, correct accuracy checking as a result, judge that whether correct the result meets the demands, if do not meet the demands, then adjusts the selection result at reference mark.This check result can be used as the geometric accuracy evaluation reference index of correcting the result.

Following portions is described the concrete technology contents according to correcting method of the present invention.

The reference mark image database makes up

Cycle and precision that the fast development of modern Remote Sensing Technical, remote sensing images are obtained also are improved gradually, and this foundation for the reference mark image database provides reliable data base.And the development of Computer Automatic Recognition technology, then the efficient use for the reference mark database provides technical foundation.According to embodiments of the invention, set up the reference mark image database, the attribute information and the image information at reference mark are unified to build library management, realized the target of " once build the storehouse, partial update is repeatedly used ".Simultaneously, the automatic matching technique of image is incorporated in the choosing automatically of reference mark, realizes that the automatic or semi-automatic geometric correction of remote sensing image is handled.

In one embodiment, each reference mark image sheet all comprises two kinds of data: view data and attribute data.Wherein attribute data is used for describing geographical location information, comprises four aspect contents:

1. the GCP geographical location information is described, as three-dimensional coordinate X, Y, Z;

2. some necessary supplementarys of geographic coordinate are described, as the coordinate system that adopted, projection pattern, ellipsoidal parameter etc.;

3. the supplementary of description control point image is as type of sensor, wave band, figure image width height, image resolution ratio etc.;

4. the feature description chosen of GCP, as the crossing or the bridge central point of road, these information can be used as the subsidiary conditions of inquiry.

The reference mark image database adopts to a large amount of reference mark image sheet that the mode of database is stored, management and service.Traditional method is the attribute at reference mark and view data to be isolated come on data storage management, what store in database is the file pointer of a correspondence image, view data is then stored separately with file mode in the database outside, this mode has been destroyed the integrality and the safeness of Data Bank of control point information, very easily loses the image information at reference mark owing to the mistake deletion of file.Adopt the data storage form that image and attribute are combined according to reference mark of the present invention image database, promptly reference mark image sheet corresponding image information is adopted binary large object BLOB (Binary Large Object) type and carries out integrated storage administration as a field of list structure with attributes such as controlling coordinate, ellipsoid type, projection pattern.

DOM(digital orthophoto map to existing 1:10000 yardstick, Digital Orthophoto Map), the DEM(digital elevation model, Digital Elevation Mode) result map is put in order, it is strong to choose the trend of the times, clean mark, the evident characteristic area, as intersection, bridge, the image sheet is gathered by a certain size (as 512 * 512 pixels) in area such as ridge angle point, from DOM, obtain planimetric coordinates information, ellipsoid, information such as projection, from the dem data of correspondence, obtain the height value of respective regions simultaneously, will obtain image information then and attribute information is unified into library storage.

Fig. 2 a is the reference mark image sheet (512 * 512 pixels, resolution are 1 meter) that the DOM image from the somewhere cuts out; Fig. 2 b be areal No. three, resource face panchromatic image (resolution is 2.1 meters), Fig. 2 c is that the resource of areal is faced multispectral image (resolution is 5.8 meters) for No. three.Contrast Fig. 2 b and Fig. 2 c can choose corresponding same place very intuitively as the reference mark from Fig. 2 a.

In addition, the coverage of reference mark image database and memory space are after certain scale, if physical sequential is routinely retrieved one by one, then the plenty of time can be expended, be unfavorable for practical application, in order to retrieve control strip fast, in actual applications, need carry out the subregion storage by geographic coordinate to the reference mark database.In the subregion process, according to the complexity of topographic relief, atural object and concrete application, to interested area and the closeer area of distribution, reference mark, reduce the division scope, and enlarge the division scope for secondary sites or water field of big area etc., the final interior number of control points in each zone that guarantees is consistent substantially, accelerates the speed of retrieval.

The retrieval of reference mark image database

Search function is to weigh the important technology index of Database Systems.For the reference mark image database, when carrying out geometric correction, can to retrieve available reference mark image fast be one of application system success or failure key of building according to waiting to correct image.

Fig. 3 is a kind of indicative flowchart of reference mark image sheet retrieval among the step S3.At step S301,, carry out the retrieval of based target regional center point longitude and latitude based on the geographical position range of the summary of being estimated among the step S1; At step S303, screen according to the attribute information of reference mark image sheet; At step S305, carry out content-based advanced search; At step S307, obtain required reference mark image sheet.Introduce the step of reference mark image sheet retrieval below respectively.

At first, according to waiting that the orbit parameter of correcting image predicts the geographic range of its summary, and then estimate the general geographic location scope of target area inner control point.Because the general location of prediction has error, there is geometry deformation in remote sensing image to be corrected in addition, generally all can be when the scope of the geographic position of estimation target area with reference to the error radius R of an estimation, its value is generally 2 to 3 times of reference mark image size.Wherein, the upper left angle point latitude and longitude coordinates of supposing this scope for (L1, B1), bottom right angle point latitude and longitude coordinates be (L2, B2), the coordinate of target area image center to be checked be (L0, B0); L1=L0-R then, L2=L0+R, B1=B0-R, B2=B0+R.

The retrieval of the center longitude of based target area image promptly is based on the retrieval in the summary geographic position of being estimated.Data item L and B in the image database record sheet of reference mark are for the reference mark video search based on general location provides pacing items.By above-mentioned coordinate relation as can be known, the expression formula for search that can set up " L1≤Li﹠﹠Li≤L2﹠﹠B2≤Bi﹠﹠Bi≤B1 " carries out relation retrieve.

Based target image attribute retrieval promptly according to the resolution of target image, sensor type, imaging time etc., filters out available reference mark image sheet, effectively dwindles the result for retrieval scope.

According to embodiments of the invention, at first estimation waits to correct the geographical position range of summary of remote sensing image target area, the scope of dwindling retrieval by the retrieval of based target regional center point longitude and latitude; Judge then available reference mark image sheet resolving range, the time reach sensor type etc., the scope of further dwindling result for retrieval mutually.Concrete retrieval mode has: based on the retrieval of given coordinate range, based on sensor, the time reach attribute information integrated retrieval such as resolution, the retrieval technique of content-based (demand distribution feature, color characteristic, shape facility, textural characteristics) mutually.To after the retrieval of reference mark image sheet, the reference mark image quantity that satisfies condition can significantly reduce, and can meet the demands basically according to this retrieval flow.

The coupling that the reference mark is right

The step of most critical is choosing of reference mark of the same name in the correction procedure of remote sensing image, and this is the key of decision image rectification robotization.And mate automatically according to existing reference mark image data is the automatic correction basis of realizing image.

An advantage of reference mark image sheet is that existing geography information has image texture information again, therefore can adopt the automatic matching algorithm of image to find the image point position of waiting to correct image of reference mark image sheet correspondence, it is right that the geographic coordinate information of taking out reference mark image sheet then just can be formed a reference mark that has comprised object coordinates and picture side's coordinate.

According to remote sensing images geometric correction method of the present invention, at reference mark image sheet with wait to correct image the time mutually, all inconsistent characteristics in aspect such as resolution, imaging angle, adopted how much invariant features to extract and matching algorithm based on SIFT, and adopt method such as coarse fuzzy C-Mean Method that the mistake match point is rejected automatically, adopt classical least square matching algorithm to carry out the coupling that becomes more meticulous of sub-pixel-level at last.

Idiographic flow as shown in Figure 4, at step S501, according to the coordinate information at reference mark, wait to correct the initial coordinate that the metadata information of image and imaging model etc. calculate corresponding picture point, then by the size of reference mark image sheet from waiting that correcting image cuts out image blocks to be searched; At step S503, utilize the Sift algorithm to reference mark image sheet and the cutting after image blocks to be searched mate, obtain preliminary matching result information; At step S505, adopt methods such as coarse fuzzy C-Mean Method and geometrical constraint to miss the rejecting of match point, it is right to keep reliably accurate match point; At step S507, utilize least-squares algorithm that matching result is carried out the essence coupling, obtain the matching precision of sub-pixel; At step S509, the reference mark that the match is successful outputed in accordance with regulations form comprise in control period, object coordinates, the control point information file as square coordinate.

As shown in Figure 1, carry out above-mentioned matching operation for a plurality of reference mark image sheet of waiting to correct in the image coverage, it is right to produce a plurality of reference mark.If the reference mark is not enough to quantity, then obtain to wait to correct more control point image sheet in the image capturing range in covering by continuing retrieval, thereby density or quantity to reference mark image sheet increase, it is right that automatic then mode of mating obtains more control point, correct when requiring when the right quantity in reference mark and distribution meeting geometric, just can carry out the geometric correction processing of robotization.Provide the relevant concrete technology contents of reference mark image sheet coupling below according to the embodiment of the invention.

The SIFT algorithm

The thought that the SIFT algorithm is selected based on the characteristics of image yardstick, set up the multiscale space of image, under different scale, detect same unique point, determine its place yardstick when determining characteristic point position, to reach the nonshrink purpose of putting of yardstick, in addition, this algorithm is rejected lower point of some contrasts and skirt response point, and extracts the invariable rotary feature descriptor to reach the purpose of anti-affined transformation.This algorithm mainly comprises: (1) sets up metric space, seeks candidate point; (2) accurately determine the key point position, reject point of instability; (3) determine the mould and the direction of key point gradient; (4) extract feature descriptor.

1. the generation of metric space

The metric space theory its objective is the multiple dimensioned feature of simulated image data when coming across computer vision field the earliest.Koendetink proves that in the literature Gaussian convolution nuclear is unique transformation kernel of realizing change of scale, and people such as Lindeberg prove further that then gaussian kernel is unique linear kernel.

Two-dimensional Gaussian function is defined as follows:

$G(x,y,\mathrm{\&sigma;})=\frac{1}{{2\mathrm{\&pi;\&sigma;}}^2}{e}^{-({\mathrm{<figure-callout\; id="2"\; label="x"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">x</figure-callout>}}^2+y^2)}/{2\mathrm{\&sigma;}}^2$

σ has represented the variance of Gauss normal distribution.

One width of cloth two dimensional image, the metric space under different scale are represented and can be obtained by image and gaussian kernel convolution:

L(x,y,σ)=G(x,y,σ)*I(x,y)

Fig. 5 shows one group of example of the bidimensional image under the different scale space.

In order effectively to detect stable key point, difference of Gaussian metric space (DoG scale-space) has been proposed at metric space.Utilize the Gaussian difference pyrene and the image convolution of different scale to generate:

D(x,y,σ)=(G(x,y,kσ)-G(x,y,σ))*I(x,y)=L(x,y,kσ)-L(x,y,σ)

The DoG operator calculates simple, is the approximate of the normalized LoG operator of yardstick.

Make up image pyramid then, image pyramid is the O group altogether, and every group has the S layer, and the image of next group is obtained as down-sampled by last set of diagrams.Fig. 6 shows one group of example of difference of Gaussian metric space (DoG) image.In Fig. 6, to metric space octave, raw video produces the image of the metric space of a series of settings through repeatedly Gaussian convolution computing.DoG image on the right is to carry out calculus of differences by the image behind the gaussian filtering that closes on to produce.After each rank, Gauss's image do the factor be 2 down-sampled, and repeat this process.

2. spatial extrema point detects

In order to seek the extreme point of metric space, each sampled point will with its all consecutive point relatively, the image area and consecutive point of scale domain than it are big or little to see it.As shown in Figure 7, middle check point and it with 9 * 2 points of 8 consecutive point of yardstick and neighbouring yardstick correspondence totally 26 points relatively, to guarantee all to detect extreme point at metric space and two dimensional image space.

3. the key point position is determined to reject with point of instability

A) the key point exact position is determined

(Taylor's secondary expansion σ) is carried out least square fitting for x, y, further determines the exact position and the yardstick of key point by the extreme value of calculating fitting surface to utilize metric space function D.Coordinate and yardstick that key point is final can be as accurate as sub-pixel-level.

With Taylor's formula launch D (σ), then the sampled point initial point is for x, y:

$D\left(X\right)=D+\frac{{\&PartialD;D}^T}{\&PartialD;X}X+\frac{1}{2}{X}^T\frac{{\&PartialD;}^{2}D}{{\&PartialD;X}^2}X$ (wherein Χ=(x, y, σ) ^T)

To the X differentiate, and to make it be zero, that is:

, just can try to achieve the sampling initial point the position

For: $\hat{X}=-\frac{{\&PartialD;}^2{D}^{-1}}{{\&PartialD;\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">X</figure-callout>}}^2}\frac{\&PartialD;D}{\&PartialD;X}$

Be $\left[\begin{array}{ccc}\frac{{\&PartialD;}^{2}D}{{\&PartialD;\mathrm{\&sigma;}}^2}& \frac{{\&PartialD;}^{2}D}{\&PartialD;\mathrm{\&sigma;y}}& \frac{{\&PartialD;}^{2}D}{\&PartialD;\mathrm{\&sigma;x}}\\ \frac{{\&PartialD;}^{2}D}{\&PartialD;\mathrm{\&sigma;y}}& \frac{{\&PartialD;}^{2}D}{{\&PartialD;\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">y</figure-callout>}}^2}& \frac{{\&PartialD;}^{2}D}{\&PartialD;\mathrm{yx}}\\ \frac{{\&PartialD;}^{2}D}{\&PartialD;\mathrm{\&sigma;x}}& \frac{{\&PartialD;}^{2}D}{\&PartialD;\mathrm{yx}}& \frac{{\&PartialD;}^{2}D}{{\&PartialD;x}^2}\end{array}\right]\left[\begin{array}{c}\mathrm{\&sigma;}\\ y\\ x\end{array}\right]=-\left[\begin{array}{c}\frac{\&PartialD;D}{\&PartialD;\mathrm{\&sigma;}}\\ \frac{\&PartialD;D}{\&PartialD;y}\\ \frac{\&PartialD;D}{\&PartialD;x}\end{array}\right]$

B) low contrast is rejected

By

If, | D (X) |＜0.03, then this contrast is lower, rejects.

C) removal of skirt response

An extreme value that defines bad difference of Gaussian operator has bigger principal curvatures in the place across the edge, and in the direction of vertical edge less principal curvatures is arranged.Principal curvatures is obtained by the Hessian matrix H of a 2x2:

$H=\left[\begin{array}{cc}{D}_{\mathrm{xx}}& D_{\mathrm{xy}}\\ D_{\mathrm{xy}}& D_{\mathrm{yy}}\end{array}\right]$

Derivative is estimated to obtain by the adjacent difference of sampled point.

The principal curvatures of D and the eigenwert of H are directly proportional, and make that α is an eigenvalue of maximum, and β is minimum eigenwert, then:

Tr(H)=D _xx+D _yy=α+β

Det(H)=D _xxD _yy-(D _xy) ²=αβ

Make α=γ β, then:

$\frac{{\mathrm{Tr}\left(H\right)}^2}{\mathrm{Det}\left(H\right)}=\frac{{(\mathrm{\&alpha;}+\mathrm{\&beta;})}^2}{\mathrm{\&alpha;\&beta;}}=\frac{{(\mathrm{r\&beta;}+\mathrm{\&beta;})}^2}{{\mathrm{r\&beta;}}^2}=\frac{{(r+1)}^2}{r}$

Value minimum when two eigenwerts equate, increase along with the increase of r.Therefore, in order to detect principal curvatures whether under certain thresholding r, only need to detect:

$\frac{{\mathrm{Tr}\left(H\right)}^2}{\mathrm{Det}\left(H\right)}<\frac{{(r+1)}^2}{r}$

Get r=10.

4. key point gradient-norm and direction calculating

Utilize the gradient direction distribution character of key point neighborhood territory pixel to be each key point assigned direction parameter, make operator possess rotational invariance.

$m(x,y)=\sqrt{{(L(x+1,y)-L(x-1,y))}^2+{(L(x,y+1)-L(x,y-1))}^2}$

$\mathrm{\&theta;}(x,y)=\mathrm{\&alpha;}\tan 2{\left((L(x,y+1)-L(x,y-1)\right)}^2+{(L(x+1,y)-L(x-1,y)))}^2$

Following formula is that (x y) locates the mould value and the direction formula of gradient.Wherein the used yardstick of L is each key point yardstick at place separately.

5. feature descriptor generates

Fig. 8 shows the process by key point neighborhood gradient information generating feature vector.The central point of Fig. 8 left half is the position of current key point.At first, coordinate axis is rotated to be the direction of key point, to guarantee rotational invariance.Next be that 8 * 8 window is got at the center with the key point.In Fig. 8, each little lattice is represented a pixel of key point neighborhood place metric space, the direction of arrow is represented the gradient direction of this pixel, and arrow length is represented the gradient-norm value, and circle is represented the scope (the pixel gradient directional information contribution the closer to key point is big more) of Gauss's weighting among the figure.On per 4 * 4 fritter, calculate the gradient orientation histogram of 8 directions then, draw the accumulated value of each gradient direction, can form a seed points, shown in Fig. 8 right half.Among this figure key point by 2 * 2 totally 4 seed points form, each seed points has 8 direction vector information.The thought of this neighborhood directivity information associating has strengthened the antimierophonic ability of algorithm, also provides fault-tolerance preferably for the characteristic matching that contains positioning error simultaneously.

In computation process, in order to strengthen the robustness of coupling, alternatively, to each key point use 4 * 4 totally 16 seed points describe, just can produce 128 data for a key point like this, promptly finally form the 128 SIFT proper vectors of tieing up.The influence that this moment, the SIFT proper vector was removed geometry deformation factors such as dimensional variation, rotation continues the length normalization method with proper vector again, then can further remove the influence of illumination variation.

Fig. 9 shows the one group example of different images by key point neighborhood gradient information generating feature vector plot.

Remote Sensing Images Matching based on the SIFT descriptor

After independent piece image is carried out above-mentioned feature extraction and feature description, just obtained feature and descriptor thereof all among this figure, being made as image 1(is realtime graphic), its unique point quantity is m.Accomplish two width of cloth images match, just obtain the pixel (referring to unique point here) that mates in two width of cloth images.At first, carry out identical feature extraction and feature description process, obtain the feature that quantity is n other piece image (image 2, i.e. reference picture); Secondly exactly will be in image 1 search out the unique point of r in n unique point in m unique point and the image 2 to mating, r≤m wherein, r≤n, and in order to calculate the geometric relationship between image 1 and the image 2 robust, should guarantee r 〉=8,, just need to adjust the precision of coupling if can not search out the r that satisfies condition, such as the similarity requirement that reduces between the matching characteristic point, feasible to obtain more matching characteristic point right.

The coupling of unique point comes down to the coupling of its descriptor, and the descriptor of unique point is that this feature has been carried out a quantitative description in fact, can be applied to matching algorithm.The coupling of unique point descriptor is carried out in describing the space in fact, be the vector of 128 dimensions such as SIFT feature description, so the SIFT descriptor mates in 128 dimension spaces.In describing the space, the matching degree of unique point descriptor is then measured with distance, and two nearest descriptors have generally just been represented a pair of unique point of coupling.And in describing the space, the definition of following two kinds of distances is arranged generally:

A) Euclidean (Euclidean) distance.Promptly be the p dimensional vector for 2 x in the p dimension space and y(), have their Euclidean distance to be defined as:

$d_E(x,y)=\sqrt{{(x_1-y_1)}^2+...+{(x_p-y_p)}^2}=\sqrt{{(x-y)}^T(x-y)}$

X=(x in the formula ₁..., x _p) ^T, y=(y ₁..., y _p) ^TEuclidean distance can be applicable to based on histogrammic descriptor, and promptly each dimension of descriptor all has identical weight, such as the SIFT descriptor, and GLOH descriptor and PCA-SIFT descriptor etc.

B) Ma Shi (Mahalanobis) distance.If each dimension of descriptor all has different weights, then need to use mahalanobis distance to measure distance between them, the weight vectors of establishing descriptor at this moment is s=(s ₁..., s _p) ^T, then 2 x in the P dimension space and the mahalanobis distance between the y are:

$d_M(x,y)=\sqrt{{\left(\frac{x_1-y_1}{s_1}\right)}^2+...+{\left(\frac{x_p-y_p}{s_p}\right)}^2}=\sqrt{{(x-y)}^T{D}^{-1}(x-y)}$

In the formula

Mahalanobis distance can be applied to the controlled filtering descriptor of direction, the constant descriptor of difference, constant descriptor of square and complex filter descriptor etc.

In an embodiment according to the present invention, use SIFT feature description to mate, it is with the standard of Euclidean distance as judge descriptor matching degree.Figure 10 shows according to an embodiment of the invention the matching algorithm flow process based on the Sift descriptor.Part among Figure 10 in the square frame is represented the extraction and the feature description process of SIFT unique point.

According to algorithm shown in Figure 10, before the output matching result, also to select match point by " minimum distance and second closely compares back-and-forth method ".Particularly, at first set a rational threshold value t, the distance of relatively being had a few in A and n the point then, find out then with A apart from the some B of minimum and with A apart from the second little some C, and as if d _ABd _ACDuring＜t, think that just A and B are reasonable correct match points.The benefit of doing like this is when A has a lot of similar coupling, promptly differs when not being a lot of with the distance of A minimum with the little distance of A second in n point, can think that this is not one and reasonably mates.Have only when the minimum distance and the second little distance are mutually far short of what is expected, when promptly the minimum point of distance " holds a safe lead " other, just admit it, because this has guaranteed that substantially this is a highly stable correct coupling rather than an ambiguous coupling.

The mistake match point is rejected

The matching characteristic point that has obtained realtime graphic and reference picture to after, basically realized the task of images match, but the unique point of coupling can only representative image between local relation, the matching characteristic point of limited quantity can not reflect fully image between holotopy.Match point centering in correspondence, may exist the bigger point of mistake coupling or matching error right, because the right existence of these points will influence matching precision, for the application based on the match point file impacts, so after initial matching finished, it was also extremely important to explore the right elimination method of mistake match point.Further, the automatic rejecting of mistake match point effectively improves the accuracy rate of mating in robotization ground.

Figure 11 shows the matching algorithm flow process that has comprised mistake match point rejecting step according to an embodiment of the invention.In Figure 11, at first carry out Feature Points Matching based on how much invariant features according to aforementioned SIFT algorithm, screen deletion for matching result then, wherein relate to coarse fuzzy C-Mean Method, random sampling consistency detecting method (RANSAC method), and the least square method of fitting of a polynomial.These method steps are described respectively below.

As previously mentioned, utilize the SIFT algorithm extract respectively real-time imaging with reference to the proper vector information of image.SIFT proper vector according to two width of cloth images, the Euclidean distance of employing key point proper vector is used as the similarity determination tolerance of key point in two width of cloth images, the minimum similarity determination tolerance as key point in two width of cloth images of the ratio (NN/SCN) of the distance that concrete principle is arest neighbors (NN) and second neighbour (SCN).Afterwards, utilize coarse fuzzy C-Mean Method to miss the preliminary rejecting of match point.

Coarse fuzzy C-Mean Method that the mistake match point is rejected

The basic thought of cluster analysis is very simple, directly perceived and simple, and it is to classify according to each pattern feature similarity degree to be classified, and the similar class that is classified as is dissimilar as an other class.Cluster analysis comprises two substances: the tolerance of pattern similarity and clustering algorithm.Fuzzy C-average (Fuzzy C-Means) algorithm is a kind of with square error in the infima species be clustering criteria, calculate the degree of membership that each sample belongs to each fuzzy subset (cluster), come implementation algorithm by the Pickard iteration between the necessary condition of objective function minimization.

Definition 1: establish X={x ₁, x ₂, x ₃..., x _nFor treating the object of classification collection, for i class w _i, its barycenter is v _i, definition

For belonging to w certainly _iThe object set of class, then:

A) if , then for $\&ForAll;k\&Element;\{\mathrm{1,2},...,n\},k\&NotEqual;j,x_j\&NotElement;{\stackrel{\&OverBar;}{R}w}_i;$

B) simultaneously in order to guarantee A _iBe unlikely to obtain too small, need to satisfy for , then exist at least k ∈ 1,2 ..., n} makes

A wherein _iBecome and go up approximate limit, go up approximate limit and portrayed the border that institute might belong to the object of i class, limit the scope that is defined if certain object does not belong to be similar to, then it belongs to the negative territory of this class, does not promptly belong to this class fully.

Definition 2: coarse fuzzy C-average (Rough Fuzzy C-Means, RFCM) objective function of algorithm is:

$J_m(U,V)=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}}_{x_j\&Element;\stackrel{\&OverBar;}{R}{w}_i}{u}_{\mathrm{ij}}^m{d}_{\mathrm{<figure-callout\; id="2"\; label="ij"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">ij</figure-callout>}}^2$

Constraint condition is:

●u _ij∈[0,1]

● $0<\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}<N$

● ${\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}}_{x_j\&Element;\stackrel{\&OverBar;}{R}{w}_i}{u}_{\mathrm{ij}}^m=1$

Utilize lagrange's method of multipliers, can obtain unconfined criterion function:

$F=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m_{x_j\&Element;\stackrel{\&OverBar;}{R}{w}_i}{d}_{\mathrm{ij}}^2-\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_i\left(\underset{i=1}{\overset{c}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}-1_{x_j\&Element;\stackrel{\&OverBar;}{R}w}\right)$

The extremum conditions of following formula is:

$\frac{\&PartialD;F}{{\&PartialD;u}_{\mathrm{ij}}}=0,\frac{\&PartialD;F}{{\&PartialD;\mathrm{\&lambda;}}_j}=0$

Following formula is calculated and can get:

$u=1/\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}{(d_{\mathrm{<figure-callout\; id="2"\; label="ij"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">ij</figure-callout>}}^2/d_{\mathrm{kj}}^2)}^{\frac{1}{m-1}}_{x_j\&Element;\stackrel{\&OverBar;}{R}{w}_i}$

The constant formula of centroid calculation formula is:

$v_i=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m{x}_j/\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{u}_{\mathrm{ij}}^m$

Can obtain following two character from the RFCM algorithm:

● $x_j\&Element;\stackrel{\&OverBar;}{R}{w}_i\&DoubleLeftRightArrow;u_{\mathrm{ij}}=1;$

● $u_{\mathrm{ij}}^{\mathrm{FCM}}\&le;u_{\mathrm{ij}}^{\mathrm{RFCM}}$

The main thought of RFCM algorithm is that an object that belongs to certain class has been divided into three set sure, possible and that negate, with square error in the infima species of all possible object be that clustering criteria carries out cluster.Different being of RFCM algorithm and FCM algorithm maximum, it thinks x _jBelong to w _iDegree of membership u _IjCalculating only comprise x in approximate with last _jClass relevant, if certain class w _iGo up and not comprise x in approximate _j, then this class is to x _jDegree of membership without any the contribution.

The objective function of RFCM algorithm be seek in each cluster infima species square distance and, it has got rid of the object that can not belong to class beyond the suprasphere.For u _IjCalculating, if for class w ₁, object x _jIt is approximate only to be positioned at going up of class , in, then

, i.e. u _1j=1, if x _jBelong to With

Common factor, u then _1j, u _2jCalculating only with class w ₁, and w ₂Relevant, and and w ₃Irrelevant, that is:

$u_{1j}=1/{(1+\frac{{\left|\right|x_j-{\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="HDA00003063937700061.png,HDA00003063937700092.png"\; state="\{\{state\}\}">v</figure-callout>}}_1\left|\right|}^2}{{\left|\right|x_j-{\mathrm{<figure-callout\; id="2"\; label="v"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">v</figure-callout>}}_2\left|\right|}^2})}^{\frac{1}{m-1}},$ $u_{1j}=1/{(1+\frac{{\left|\right|x_j-{\mathrm{<figure-callout\; id="2"\; label="v"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">v</figure-callout>}}_2\left|\right|}^2}{{\left|\right|x_j-{\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="HDA00003063937700061.png,HDA00003063937700092.png"\; state="\{\{state\}\}">v</figure-callout>}}_1\left|\right|}^2})}^{\frac{1}{m-1}}$

If x _jDo not exist

In, u then _1j=0.

Objective function

, provided u _IjMore new formula

According to this formula, with v _jFor center radius the u in the suprasphere of α _IjBe 1.

Based on above-mentioned mathematical analysis, reject step and mainly comprise according to the coarse fuzzy C under the geometrical constraint of the embodiment of the invention-Mean Method mistake match point:

To reaching the resolution information of image, calculate the geometrical constraint condition between image according to the match point that generates:

● under the resolution situation, ask the right slope k of corresponding point and apart from S;

● under the different resolution situation, ask corresponding point to the intersecting point coordinate of line extended line (X, Y).

Then, utilize coarse fuzzy C-Mean Method (RFCM), ask (k, S) or (to carrying out cluster analysis, it is right to delete inhomogeneous point to all match points for X, degree of membership Y), only keeps the point that is included in the sure gathering set.

I. determine that class counts c (2≤c≤N), parameter m, initial matrix, the last approximate boundaries A of class _iWith a suitable decimal, s=0;

Ii. calculate barycenter

Iii. if

U then _Ij=0, otherwise upgrade

Iv. if || U ^s-U ^S+1||＜ε then stops, otherwise s=s+1 forwards b to.

As shown in figure 11, set number (for example, 7) match point then adopts the random sampling consistency detecting method to carry out further match point and screens out if carry out still having after the match point examination surpassing through coarse fuzzy C-Mean Method.

Random sampling consistency detecting method (RANSAC method)

A) calculated amount of RANSAC algorithm

In the RANSAC algorithm, require to guarantee under certain fiducial probability that having the data of one group of sampling at least in the sampling of M group is interior point (inliers) entirely.Utilize the following formula can be in the hope of the minimum sampling number M that meets the demands.

1-(1-(1-ε) ^m) ^M=p

Wherein: ε is data error rate (exterior point (outliers) is in the shared ratio of raw data), the minimum data amount that m needs for the computation model parameter, and P is a fiducial probability.As can be seen from the above equation, M and ε, m, P is exponential relationship.When following table showed P=0.95, M was with m, the situation that ε changes.

As can be seen from the above table, when the model more complicated, when ε is higher, M is very big, directly causes the RANSAC efficiency of algorithm to descend.Need the time if from raw data, randomly draw one group of sampling; Need time t by one group of sample calculation model parameter; With the data testing model parameter average out to t that takes time, then need time N with N data checks (all data check) _tTherefore the RANSAC algorithm is required computing time:

T=M(t _s+t _c)+MN _t

Wherein, M (t _s+ t _c) be that sampling extraction of M group and model parameter are calculated the time that needs, MN _tCheck the time that needs for M model parameter.

Find out that from this formula the time that the RANSAC algorithm needs is made up of two parts:

● the time that the sampling of M group is selected and model parameter estimation needs;

● the time that the model parameter check needs.Under the situation that model is determined, data error rate is determined,

In order to guarantee result's fiducial probability, M can not reduce.

Therefore, in order to improve the efficient of algorithm, can only be from reducing the model parameter quantity that participates in check, the time that the check of minimizing model parameter needs.

B) RANSAC algorithm steps

Calculate minimum sampling number M according to fiducial probability P and data error rate ε; Calculate the corresponding model parameter of sampling,, obtain the inliers quantity of each model parameter with all raw data testing model parameter qualities; Select the optimization model parameter according to the variance of inliers quantity and error; Inliers with optimization model parameter correspondence estimates the final mask parameter.

The smart coupling of least square method

Carry out after further match point is selected through the RANSAC algorithm, can also utilize the least square method of fitting of a polynomial reject the match residual error greater than

The match point of error is right doubly.

As shown in figure 11, less if coupling is counted, for example be not more than 7, alternatively, can adopt parameter adaptive SIFT algorithm to increase for the match point that screens.

Satellite remote-sensing image is owing to extraneous factors such as being subjected to weather, sunlight, block has a strong impact on, and there is difference because of the imaging pattern camera parameter model trajectory between problem such as image translation that factors such as different imaging times, angle, distance cause, rotation, convergent-divergent and the various sensor, even adopt the SIFT algorithm, few or the basic extraction of the unique point that also may occur extracting under extreme case is less than the phenomenon of unique point, causes the mistake coupling or it fails to match.At this situation, can adopt improved parameter adaptive SIFT image matching method, according to different imaging characteristics and quality, choose corresponding strategy and determine the respective threshold that unique point is surveyed according to relevant weight definition rule.

The parameter that relates in SIFT feature extraction and matching process is more, for the effect of verifying that each parameter plays in feature extraction and matching process, carried out the coupling experiment, main by revising the size of certain parameter to be verified, promptly this parameter is regulated according to theoretical rule change in the algorithm, meanwhile keep other parameter constants, the match point number during by each parameter size of coupling counting statistics, mistake match point number, match time etc., and carry out statistical study.

By to the experiment the statistics analysis as can be known, its general rule is as follows:

● when the sampling interval of every rank (Octave) increased gradually, the match point number increased rapidly, and mistake match point quantity does not change, and increase match time;

● increase Gaussian convolution nuclear σ gradually, the coupling whole downtrending that becomes of counting, but the result is best when σ=1.6;

● contrast threshold increases gradually, and the match point number reduces gradually, and mistake match point data reduce gradually;

● along with the increase gradually of bent ratio threshold value, the match point number increases gradually, and the mistake coupling is counted and become hardly;

● coupling threshold value minimum distance increases with the ratio of time minimum distance, and the match point number increases gradually, and the mistake coupling is counted and also progressively increased.

Based on above-mentioned experimental result, parameter adaptive SIFT algorithm has been proposed.The step of parameter adaptive SIFT algorithm is to determine the threshold value of each parameter automatically by the information calculations of image self, thereby by changing the unique point quantity that parameter plays increases extraction, the purpose of increase match point quantity.Specifically can comprise the steps:

A) calculate the average gray of input image or image done Auto Laves, if then with average gray divided by 10 results less than 1.0, then reduce the threshold value of contrast, increase the unique point quantity of extracting;

B) satisfying under the situation of matching precision,, then increasing bent ratio threshold value if coupling was counted less than 7 o'clock;

C) if through a) b) the still less or number of matches of the feature quantity extracted of two steps more after a little while, under the situation that does not influence matching precision, suitably amplify the proportion threshold value of minimum distance and time minimum distance.

As shown in Figure 1, after having carried out the reference mark coupling, according to sensor type and imaging characteristics, select corresponding geometric correction model (step 9), then based on the reference mark of being mated to and selected geometric correction modelling pixel coordinate and terrestrial coordinate between transformational relation (step 11) is carried out geometric correction (step 13) based on the transformational relation between digital elevation model (DEM) and pixel coordinate and the terrestrial coordinate to raw video at last.Below the main selection of geometric correction model and the step of geometric correction described.

Geometric correction model and selection thereof

Linear array push-broom type imaging sensor be line by line with the time sequential mode obtain two dimensional image.Generally be to form a bar chart picture earlier on image planes, satellite is pushed ahead along predetermined track then, and the scanning back forms a width of cloth bidimensional image one by one, and imaging mode as shown in figure 12.Each row pixel is in the synchronization imaging and for central projection on the image, and whole image is the multicenter projection.

In Figure 12: p _kBe arbitrary picture point on the image, x _kBe the x coordinate of the last imaging point of sweep trace k, c is a sensor master distance, O _kBe the projection centre of sweep trace k, o _kBe the principal point of sweep trace k, l _kFor sweep trace k from projection centre O _kThe light that sends.

Problem the most basic will be set up rational remote sensing image imaging model exactly in the geometric correction, and the imaging model of so-called remote sensing image is meant the coordinate set up on the image (x, y) topocentric terrestrial coordinate (X, Y, Z) mathematical relation between corresponding with it.That is:

X=f _x(x,y,g)

Y=f _y(x,y,h)

G wherein, h is the influence of other factors.

The imaging model of remote sensing image can be divided into two big classes: physical model and universal model.

Physical model be meant be considered to as the time cause factors such as the physical significance of deformation of image such as surface relief, earth curvature, atmospheric refraction, camera distortion, utilize these physical conditions to be construed as the picture geometric model then, the most representational is that collinearity condition equation is the strict imaging model on basis.

Universal model is not considered imaging mechanism, but directly describes the geometric relationship of picture point and object point with mathematical function, and it has the characteristics of generality, confidentiality, high efficiency.Universal imaging model commonly used has polynomial expression, direct linear transformation, affined transformation, rational function model etc.;

In an embodiment according to the present invention, with strict imaging model, rational function model and multinomial model be the optional model of geometric correction.According to sensor type and imaging characteristics, select corresponding geometric correction model, for example, the SPOT5 satellite image, the strict imaging model of general selection, then select rational function model for satellite images such as IKONOS, P5, RapidEye, nothing money rail parameter or rational function model coefficient file, multinomial model then selected.

Strict imaging model is described when being imaging picture point and the strict one to one geometric process of object space point are representative with the collinearity condition equation of using always.According to one embodiment of present invention, the high-precision satellite almanac data, attitude data and the time data that provide according to SPOT5, by imaging process, set up picture point (l, p) with ground point (X, Y, Z) or

Between mathematical relation.Figure 13 shows picture point and topocentric mathematical relation.

The used coordinate system of strict imaging model may comprise following coordinate system, as shown in figure 14:

(also the claiming body coordinate system) O of navigation coordinate system ₁-X ₁Y ₁Z ₁: be that wherein initial point is at satellite barycenter, X by the decision of the attitude control system (ACS) of satellite ₁Axle, Y ₁Axle, Z ₁Axle is got three principal moments axles of satellite respectively.X ₁Axle is along satellite transverse axis, Y ₁Axle points to satellite flight direction (not exclusively overlapping), Z along the longitudinal axis ₁Axle is determined according to the right-handed system rule.The attitude of satellite is measured and carried out in body coordinate system, and three parameters describing its spatial attitude are: pitching (pitch) is around body coordinate system X ₁The rotation of axle, roll (roll) is around body coordinate system Y ₁The rotation of axle, course deviation (yaw) is around body coordinate system Z ₁The rotation of axle.

(ψ _x) _p(ψ _y) _pThe inclination angle of direction along ng a path and vertical track direction when representing P components and parts imaging of charge-coupled device (CCD) respectively.

Orbital coordinate system: be the transition coordinate system that navigation coordinate system and earth axes transform, the initial point of orbital coordinate system is at satellite barycenter, Z ₂It is reverse that axle points to the earth's core, Y ₂Axle points to the direction of satellite motion, X on the satellite orbit face ₂Axle is determined according to the right-handed system rule.

O ₂-X ₂Y ₂Z ₂Each direction by satellite in t position constantly

And speed

Two direction vectors are determined, and are as follows:

${\stackrel{\&RightArrow;}{Z}}_2=\frac{\stackrel{\&RightArrow;}{P}\left(t\right)}{\left|\right|\stackrel{\&RightArrow;}{P}\left(t\right)\left|\right|}$

${\stackrel{\&RightArrow;}{X}}_2=\frac{\stackrel{\&RightArrow;}{V}\left(t\right)\mathrm{\&Lambda;}{\stackrel{\&RightArrow;}{Z}}_2}{\left|\right|\stackrel{\&RightArrow;}{V}\left(t\right)\mathrm{\&Lambda;}{\stackrel{\&RightArrow;}{Z}}_2\left|\right|}$

${\stackrel{\&RightArrow;}{Y}}_2={\stackrel{\&RightArrow;}{Z}}_2\mathrm{\&Lambda;}{\stackrel{\&RightArrow;}{X}}_2$

Earth axes: the position of SPOT satellite and speed all are to be defined in earth agreement reference frame (ITRS), ITRF is the realization of ITRS, be that earth reference frame portion by the IERS central office is responsible for setting up and safeguarding, the station coordinates of this reference frame and speed are to utilize such as very long baseline interferometry(VLBI (VLBI), Laser Measuring defends etc. that data computation that space geodesy technique gathers comes out, it is a reference frame that constantly changes, that adopt at present is ITRF90, calculate the coordinate system that generally adopts the WGS84 coordinate system to replace the ITRF framework to determine for convenience.This be because, WGS84 is the earth's core body-fixed coordinate system, ITRF then be an inertial coordinates system, final mapping needs be body-fixed coordinate system, and ITRS and WGS84 only differ about 1 meter, the conversion parameter from ITRF90 to WGS84 is as shown in the table:

T1(m)	T2(m)	T3(m)	D(ppm)	R1(arcsec.)	R2(arcsec.)	R3(arcsec.)
							0.060	-0.517	-0.223	-0.011	0.0183	-0.0003	0.0070

ITRF90 is to the conversion parameter table of WGS84

Conversion formula is as follows:

$\left(\begin{array}{c}{X}_{\mathrm{WGS}84}\\ Y_{\mathrm{WGS}84}\\ Z_{\mathrm{WGS}84}\end{array}\right)=\left(\begin{array}{c}{T}_1\\ T_2\\ {\mathrm{<figure-callout\; id="3"\; label="T"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">T</figure-callout>}}_3\end{array}\right)+(1+D)\left(\begin{array}{ccc}0& {-R}_3& R_2\\ {\mathrm{<figure-callout\; id="3"\; label="R"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">R</figure-callout>}}_3& 0& -R_1\\ -R_2& {\mathrm{<figure-callout\; id="1"\; label="R"\; filenames="HDA00003063937700061.png,HDA00003063937700092.png"\; state="\{\{state\}\}">R</figure-callout>}}_1& 0\end{array}\right)\left(\begin{array}{c}{X}_{\mathrm{ITRF}}\\ Y_{\mathrm{ITRF}}\\ Z_{\mathrm{ITRF}}\end{array}\right)$

The WGS84 coordinate system is a measurement reference system of being made up of global the earth's core coordinate frame of reference and one group of corresponding model and WGS84 geoid surface.The initial point of this coordinate system is positioned at earth centroid, and Z axle and IERS point to identical with reference to the utmost point, and X-axis is pointed to the IERS guide meridian and crossed initial point and perpendicular to the intersection point of Z axial plane, Y-axis is determined by the right-handed system rule;

The imaging process derivation (l, p) → (X, Y, Z) as follows:

ⅰ, the pairing imaging time of calculating l capable p row pixel

t=t _c+lsp×(l-l _c)

Wherein, t _cBe the photography time of central scan row, lsp is the sweep time of every row, l _cBe the row of central scan row number

ⅱ, interpolation go out the t position of satellite constantly

$\stackrel{\&RightArrow;}{P}\left(t\right)=\underset{j=1}{\overset{8}{\mathrm{\&Sigma;}}}\frac{\stackrel{\&RightArrow;}{P}\left(t_j\right)\&times;\underset{\underset{i\&NotEqual;j}{i=1}}{\overset{8}{\mathrm{\&Pi;}}}(t-t_i)}{\underset{\underset{i\&NotEqual;j}{i=1}}{\overset{8}{\mathrm{\&Pi;}}}(t-t_i)}$

$\stackrel{\&RightArrow;}{V}\left(t\right)=\underset{j=1}{\overset{8}{\mathrm{\&Sigma;}}}\frac{\stackrel{\&RightArrow;}{V}\left(t_j\right)\&times;\underset{\underset{i\&NotEqual;j}{i=1}}{\overset{8}{\mathrm{\&Pi;}}}(t-t_i)}{\underset{\underset{i\&NotEqual;j}{i=1}}{\overset{8}{\mathrm{\&Pi;}}}(t-t_i)}$

ⅲ, the imaging inclination angle (ψ by p row _x) _p(ψ _y) _pCalculate (l, p) the imaging direction vector under navigation coordinate system

${\stackrel{\&RightArrow;}{u}}_1^{\&prime;}=\left(\begin{array}{c}-\tan \left({\mathrm{\&psi;}}_Y\right)\\ +\tan \left({\mathrm{\&psi;}}_X\right)\\ -1\end{array}\right)$

${\stackrel{\&RightArrow;}{u}}_1=\frac{{\stackrel{\&RightArrow;}{u}}_1^{\&prime;}}{\left|\right|{\stackrel{\&RightArrow;}{u}}_1^{\&prime;}\left|\right|}$

Imaging direction vector under ⅳ, the calculating orbital coordinate system

${\stackrel{\&RightArrow;}{Z}}_2=\frac{\stackrel{\&RightArrow;}{P}\left(t\right)}{\left|\right|\stackrel{\&RightArrow;}{P}\left(t\right)\left|\right|}$

${\stackrel{\&RightArrow;}{X}}_2=\frac{\stackrel{\&RightArrow;}{V}\left(t\right)\mathrm{\&Lambda;}{\stackrel{\&RightArrow;}{Z}}_2}{\left|\right|\stackrel{\&RightArrow;}{V}\left(t\right)\mathrm{\&Lambda;}{\stackrel{\&RightArrow;}{Z}}_2\left|\right|}$

${\stackrel{\&RightArrow;}{Y}}_2={\stackrel{\&RightArrow;}{Z}}_2\mathrm{\&Lambda;}{\stackrel{\&RightArrow;}{X}}_2$

Wherein, With

Be respectively pixel (l, p) imaging constantly the position and the velocity interpolation of t satellite go out the t attitude angle of satellite constantly:

$a_p\left(t\right)=a_p\left(t_i\right)+(a_p\left(t_{i+1}\right)-a_p\left(t_i\right))\frac{t-t_i}{t_{i+1}-t_i}$

$a_r\left(t\right)=a_r\left(t_i\right)+(a_r\left(t_{i+1}\right)-a_r\left(t_i\right))\frac{t-t_i}{t_{i+1}-t_i}$

$a_y\left(t\right)=a_y\left(t_i\right)+(a_y\left(t_{i+1}\right)-a_y\left(t_i\right))\frac{t-t_i}{t_{i+1}-t_i}$

a _p(t _i), a _r(t _i), a _y(t _i) be respectively t _iThe attitude angle of satellite can be obtained from auxiliary data file constantly.

With the imaging direction vector under the navigation coordinate system

Be converted into the imaging direction vector under the orbital coordinate system

${\stackrel{\&RightArrow;}{u}}_2=\frac{{\stackrel{\&RightArrow;}{u}}_2^{\&prime;}}{\left|\right|{\stackrel{\&RightArrow;}{u}}_2^{\&prime;}\left|\right|}$

${\stackrel{\&RightArrow;}{u}}_2^{\&prime;}=M_p\&CenterDot;M_r\&CenterDot;M_y\&CenterDot;{\stackrel{\&RightArrow;}{u}}_1$

$M_p=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(a_p\left(t\right)\right)& \sin \left(a_p\left(t\right)\right)\\ 0& -\sin \left(a_p\left(t\right)\right)& \cos \left(a_p\left(t\right)\right)\end{array}\right]$

$M_r=\left[\begin{array}{ccc}\cos \left(a_r\left(t\right)\right)& 0& -\sin \left(a_r\left(t\right)\right)\\ 0& 1& 0\\ \sin \left(a_r\left(t\right)\right)& 0& \cos \left(a_r\left(t\right)\right)\end{array}\right]$

$M_y=\left[\begin{array}{ccc}\cos \left(a_y\left(t\right)\right)& -\sin \left(a_y\left(t\right)\right)& 0\\ \sin \left(a_y\left(t\right)\right)& \cos \left(a_y\left(t\right)\right)& 0\\ 0& 0& 1\end{array}\right]$

Imaging direction vector under ⅴ, the calculating earth axes

${\stackrel{\&RightArrow;}{u}}_3=\left[\begin{array}{ccc}{\left(X_2\right)}_X& {\left(Y_2\right)}_X& {\left(Z_2\right)}_X\\ {\left(X_2\right)}_Y& {\left(Y_2\right)}_Y& {\left(Z_2\right)}_Y\\ {\left(X_2\right)}_Z& {\left(Y_2\right)}_Z& {\left(Z_2\right)}_Z\end{array}\right]\&CenterDot;{\stackrel{\&RightArrow;}{u}}_2$

Wherein, X ₂, Y ₂, Z ₂Be orbital coordinate system O-X ₂Y ₂Z ₂Three change in coordinate axis direction vector ⅵ, topocentric location:

As shown in figure 15, measurement vector u ₃With the intersection point on ground be M, the position vector of satellite is

M connects corresponding vector with the earth's core

Because the M point is a+h at major axis, minor axis is on the ellipsoid of b+h again, so have behind the simultaneous equations:

$\stackrel{\&RightArrow;}{O_{3}M}=\stackrel{\&RightArrow;}{P}\left(t\right)+\mathrm{\&mu;}\&times;{\mathrm{<figure-callout\; id="3"\; label="u"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">u</figure-callout>}}_3\&DoubleRightArrow;\left\{\begin{array}{c}X=X_P+\mathrm{\&mu;}\&times;{\left(u_3\right)}_X\\ Y=Y_P+\mathrm{\&mu;}\&times;{\left(u_3\right)}_Y\\ Z=Z_P+\mathrm{\&mu;}\&times;{\left(u_3\right)}_Z\end{array}\right.$

$\frac{{\mathrm{<figure-callout\; id="2"\; label="X"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">X</figure-callout>}}^2+Y^2}{{\mathrm{<figure-callout\; id="2"\; label="A"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">A</figure-callout>}}^2}+\frac{Z^2}{{\mathrm{<figure-callout\; id="2"\; label="B"\; filenames="HDA00003063937700061.png"\; state="\{\{state\}\}">B</figure-callout>}}^2}=1$

Wherein,

$\left\{\begin{array}{c}A=a+h\\ B=b+h\end{array}\right.$

In conjunction with several transformation relations of front, write as the form of similar collinearity condition equation:

${\mathrm{<figure-callout\; id="1"\; label="u"\; filenames="HDA00003063937700061.png,HDA00003063937700092.png"\; state="\{\{state\}\}">u</figure-callout>}}_1=\left[\begin{array}{c}-\tan {\left({\mathrm{\&psi;}}_y\right)}_p\\ +\tan {\left({\mathrm{\&psi;}}_x\right)}_p\\ -1\end{array}\right]=\frac{1}{u}{R}_{21}^{-1}{.R}_{32}^{-1}.\left[\begin{array}{c}X-X_P\\ Y-Y_P\\ Z-Z_P\end{array}\right]$

The relation of the imaging direction vector under the acquisition navigation coordinate system and the observed ray of earth axes; Projection centre is (X _P, Y _P, Z _P), coordinate vector be (Z), projection coefficient is u for X, Y, and rotation matrix is:

$R=R_{21}^{-1}[a_r\left(t\right),a_p\left(t\right),a_y\left(t\right)]R_{32}^{-1}[P\left(t\right),V\left(t\right)]=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right]$

Can get collinearity equation:

$\frac{r_{11}(X-X_P)+r_{12}(Y-Y_p)+r_{13}(Z-Z_p)}{r_{31}(X-X_P)+r_{32}(Y-Y_p)+r_{33}(Z-Z_p)}-\tan {\left({\mathrm{\&psi;}}_y\right)}_p=0$

$\frac{r_{21}(X-X_P)+r_{22}(Y-Y_p)+r_{23}(Z-Z_p)}{r_{31}(X-X_P)+r_{32}(Y-Y_p)+r_{33}(Z-Z_p)}+\tan {\left({\mathrm{\&psi;}}_x\right)}_p=0$

Can obtain the coordinate on the image in the remote sensing image imaging process (x, y) topocentric terrestrial coordinate (X, Y, Z) the mathematics mapping relations between corresponding thus with it.

RFM model (Rational Function Model, rational function model) is a kind of in the general purpose transducer model, its essence is to describe strict imaging process by complicated mathematical model more, sets up the direct corresponding relation of object space coordinate and picture side's picpointed coordinate.The RFM model can effectively encapsulate concerning security matters parameters such as satellite orbit parameter, attribute sensors, and the user can not know under the prerequisite of relevant information, directly just can calculate the volume coordinate of picture point according to RFM model and parameter, thereby carry out various application.The RFM model generally adopts RPC parameter (Rational Polynomial Coefficients, rational polynominal coefficient) to represent that the RPC parameter generally has 80.

A, the RFM modular concept of utilizing RPC to locate

The RFM model be with picpointed coordinate (s, l) be expressed as with millet cake volume coordinate correspondingly (h) (wherein φ is a longitude, and λ is a latitude, and h is an elevation) is the polynomial ratio of independent variable for φ, λ:

$s=\frac{\mathrm{Nums}(U,V,W)}{\mathrm{Dens}(U,V,W)}$

$l=\frac{\mathrm{Numl}(U,V,W)}{\mathrm{Denl}(U,V,W)}$

Wherein,

$U=(\mathrm{\&phi;}-{\mathrm{\&phi;}}_0)/{\mathrm{\&phi;}}_s$ $V=(\mathrm{\&lambda;}-{\mathrm{\&lambda;}}_0)/{\mathrm{\&lambda;}}_s$

$W(h-h_0)/h_s$

$s=(S-S_0)/S_s$

$l=(L-L_0)/L_s$

(U, V, W) be after the standardization the ground point volume coordinate (φ, λ, h), (s, l) be after the standardization picpointed coordinate (S, L), a _i, b _i, c _i, d _i(i=1～20), φ ₀, φ _s, λ ₀, λ _s, h ₀, h _s, S ₀, S _s, L ₀, L _sBe the RPC parameter of photograph, φ ₀, λ ₀, h ₀, S ₀, L ₀Be standardization translation parameters, φ _s, λ _s, h _s, S _sL _sBe the standardization scale parameter.They are provided with image data.

Nums(U,V,W)=a ₁+a ₂V+a ₃U+a ₄W+a ₅VU+a ₆VW+a ₇UW+a ₈V ²

+a ₉U ²+a ₁₀W ²+a ₁₁VUW+a ₁₂V ³+a ₁₃VU ²+a ₁₄VW ²

+a ₁₅V ²U+a ₁₆U ³+a ₁₇UW ²+a ₁₈V ²W+a ₁₉U ²W+a ₂₀W ³

Dens (U, V, W), Numl (U, V, W), (U, V W) can be expressed as about (U, V, the form of cubic polynomial W) Denl too.

Finding the solution of b, RPC parameter

With the RFM model deformation be:

F _X=Num _s(P,L,H)-X*Den _s(P,L,H)=0

F _Y=Num _L(P,L,H)-Y*Den _L(P,L,H)=0

Then error equation is:

V=Bx-l,W

In the following formula,

$B=\left[\begin{array}{cccc}\frac{{\&PartialD;F}_X}{{\&PartialD;a}_i}& \frac{{\&PartialD;F}_X}{{\&PartialD;b}_j}& \frac{{\&PartialD;F}_X}{{\&PartialD;c}_i}& \frac{{\&PartialD;F}_X}{{\&PartialD;d}_j}\\ \frac{{\&PartialD;F}_Y}{{\&PartialD;a}_i}& \frac{{\&PartialD;F}_Y}{{\&PartialD;b}_j}& \frac{{\&PartialD;F}_Y}{{\&PartialD;c}_i}& \frac{{\&PartialD;F}_Y}{{\&PartialD;d}_j}\end{array}\right],$ (i=1…20,j=1…20)

$1=\left[\begin{array}{c}-F_X^0\\ -F_Y^0\end{array}\right]$

$x={\left[\begin{array}{cccc}{a}_i& b_j& c_i& d_j\end{array}\right]}^T$

W is a weight matrix.

According to the least square adjustment principle, error equation can be transformed to normal equation:

(B ^TB)x=B ^Tl

When the employing second order that is used to resolve the RPC model parameter or second order were above, there was overparameterization in the RPC model, in the RPC model variation of denominator very violent, so just cause design matrix (B ^TB) state variation, design matrix becomes singular matrix, causes least square adjustment not restrain.

Adopt the method for spectrum Modified Iteration Method solving method equation, no matter normal equation is good attitude, morbid state or rank defect, it resolves program and all need not add anyly when normal equation is good attitude, just can converge to exact solution through iteration several times.When the normal equation wan, speed of convergence is slow slightly, and estimated result does not have partially.Concrete formulae express is:

Be provided with (B ^TB) x=B ^TL adds x simultaneously with top formula both sides, obtains

(B ^TB+E)x=B ^Tl+x

Wherein, E is t rank unit matrixs, because unknown parameter x is all contained on the formula both sides, so can only adopt the method for iteration to find the solution, its iterative formula is:

x ^(k)=(B ^TB+E) ^-1(B ^Tl+x ^(k-1))

Multinomial model is the space geometry process of avoiding imaging, and the coordinate that directly adopts a number of times appropriate mathematical polynomial expression to describe between the image respective point of correction front and back concerns.

General polynomial is corrected transformation for mula:

u=a ₀₀+a ₁₀X+a ₀₁Y+a ₂₀X ²+a ₁₁XY+a ₀₂Y ²+a ₃₀X ³+a ₂₁X ²Y+a ₁₂XY ²+a ₀₃Y ³+...

v=b ₀₀+b ₁₀X+b ₀₁Y+b ₂₀X ²+b ₁₁XY+b ₀₂Y ²+b ₃₀X ³+b ₂₁X ²Y+b ₁₂XY ²+b ₀₃Y ³+...

In the formula, (u v) is the picture planimetric coordinates of picture point; (X Y) is its terrestrial coordinate of millet cake accordingly; a _Ij, b _IjBe polynomial coefficient.Here polynomial exponent number generally is not more than three times, because more the polynomial expression of high-order often can not improve precision and can cause the relevant of parameter on the contrary, causes the reduction of model orientation accuracy.

Based on the reference mark of coupling to the transformational relation of geometric correction modelling pixel coordinate of selecting and terrestrial coordinate after, at step S13, carry out geometric correction based on the transformational relation between digital elevation model and pixel coordinate and the terrestrial coordinate, obtain digital orthoimage through orthorectify.The elevation model here obtains from existing basic surveying and mapping result figure.

In addition, alternatively, after step S13, correct accuracy checking as a result, judge that whether correct the result meets the demands, if do not meet the demands, then carries out the manual intervention adjustment to the distribution and the quantity selection result at reference mark.

According to embodiments of the invention, when stating the satellite remote sensing images fast geometric correcting method based on the reference mark image database in realization, can adopt based on the parallel accelerated processing method of the geometric correction of OpenMP.

Handle based on parallel acceleration of the geometric correction of OpenMP

The standard of OpenMP is initiated by SGI, and it is a kind of multiprocessor multi-threaded parallel programming language towards shared drive and distributed shared memory.OpenMP has good portability, supports Fortran and C/C++ programming language, and unix system and Windows system are then supported in the operating system platform aspect.The importance of OpenMP is that it can provide a kind of simple method for writing multithread programs, need not the programmer carry out complicated threads create, synchronously, load balance and destruction work.OpenMP is for For loop statement particularly suitable, and it can be crossed by add the parallel processing that multinuclear just realized in several simple statements in original program, makes full use of the cpu resource of computer, realizes the acceleration of calculating.Therefore, it will be very effective adopting OpenMP to carry out that the unit multi-core parallel concurrent handles at the geometric correction of piecemeal, and hardware is not had too much requirement yet.

Remote sensing image is actually from a geometric space at correction procedure and is transformed into another geometric space, in case set up the correction model, then this conversion is unique, and is only relevant with the position, in addition because imaging is continuous, it also is continuous correcting the back image, therefore can by calculating four coordinates before and after the angle point conversion, set up simple relatively affined transformation model in the piece then by the thought of piecemeal, thereby simplify the calculating of correcting model, reduce calculated amount.Consider the influence of landform, the size of piecemeal is to need to consider, should be that the big local piecemeal of topographic relief is little in theory, the little local piecemeal that rises and falls is big, and such as resolution such as SPOT5, P5 all about 2.5 meters, so finally divide block size to be decided to be 15 * 15, divide block size then to be decided to be 5 * 5 for HJ_1A 1B.On this basis, in order to consider the needs to big Flame Image Process, can carry out the branch block operations to original image, be about to original image and be divided into some, each piece is regarded piece image as and corrected separately and write disk then, and is synthetic at last.

Following table shows at the different images type carries out the example handled based on the geometric correction of OpenMP.

As can be seen from the above table,, adopt piecemeal to add parallel processing and can significantly improve counting yield for strict imaging model, speed-up ratio about 50, and for rational function model and its speed-up ratio of polynomial expression also about 4.

According to the remote sensing images correcting method of the embodiment of the invention, the automatic matching technique of image is incorporated in the choosing automatically of reference mark, can realize that the automatic or semi-automatic geometric correction of remote sensing image is handled.Particularly, this method is obtained the reference mark matched data based on the reference mark image database by the automatic coupling of query and search, reference mark image for image storehouse, reference mark; Utilize the geometric correction of multiple models such as strict imaging model, rational function model, multinomial model again; Realize complete remote sensing image automatic geometric correction solution, effectively utilized existing reference mark resource, improved production efficiency.Simultaneously, this method can be calculated fast by the unit multi-core parallel concurrent, thereby further improves the correction efficient of remote sensing image.

The above only is an exemplary embodiment of the present invention, but not is used to limit protection scope of the present invention, and protection scope of the present invention is determined by appended claim.

Claims (10)

1. the satellite remote sensing images fast geometric correcting method based on the reference mark image database is characterized in that, comprising:

Determine to wait to correct the geographic range of image;

At determined geographic range of waiting to correct image, all satisfactory reference mark of retrieval in the image database of reference mark;

Automatically mate reconnaissance, the reference mark that is identified for geometric correction is right;

Judge the quantity that the reference mark of described coupling is right and distribute whether meeting geometric is corrected requirement, if, then enter next procedure, if not, described searching step then returned;

According to sensor type and imaging characteristics, select corresponding geometric correction model;

Based on the reference mark of described coupling to the geometric correction modelling pixel coordinate of described selection and the transformational relation of terrestrial coordinate;

Carry out geometric correction based on the transformational relation between digital elevation model and described pixel coordinate and the terrestrial coordinate, obtain digital orthoimage through orthorectify.

2. remote sensing images geometric correction method according to claim 1, it is characterized in that, in described geometric correction Model Selection, for the SPOT5 satellite image, select strict imaging model, select rational function model for IKONOS, P5, RapidEye satellite image,, select multinomial model for the satellite image of not supporting rail parameter or rational function model coefficient file.

3. remote sensing images geometric correction method according to claim 1 is characterized in that, described retrieval comprises:

Estimation waits to correct the geographical position range of summary of image;

Based on the geographical position range of the summary of described estimation, carry out the retrieval of based target regional center point longitude and latitude;

Attribute information according to reference mark image sheet screens;

Carry out content-based advanced search.

4. remote sensing images geometric correction method according to claim 3, it is characterized in that described attribute information according to reference mark image sheet screens the resolution, sensor type, the imaging time that comprise according to target image and screens available reference mark image sheet.

5. remote sensing images geometric correction method according to claim 3 is characterized in that, described content-based advanced search comprises based on demand distribution feature, color characteristic, shape facility, textural characteristics coming reference mark image sheet is retrieved.

6. remote sensing images geometric correction method according to claim 1 is characterized in that, described automatic coupling comprises:

According to the coordinate information at reference mark, wait to correct the metadata information of image and the initial coordinate that imaging model calculates corresponding picture point, then by the size of reference mark image sheet from waiting that correcting image cuts out image blocks to be searched;

Image blocks to be searched after utilizing the Sift algorithm to reference mark image sheet and cutting is mated, and obtains preliminary matching result information;

Adopt coarse fuzzy C-Mean Method and geometrical constraint method to miss the rejecting of match point, it is right to keep reliably accurate match point;

Utilize least-squares algorithm that matching result is carried out the essence coupling, obtain the matching precision of sub-pixel;

The reference mark that the match is successful outputed in accordance with regulations form comprise in control period, object coordinates, the control point information file as square coordinate.

7. remote sensing images geometric correction method according to claim 6 is characterized in that, the described Sift of utilization algorithm mates and comprises:

Set up metric space, seek candidate point;

Accurately determine the key point position, reject point of instability;

Determine the mould and the direction of key point gradient;

Extract feature descriptor.

8. remote sensing images geometric correction method according to claim 1 is characterized in that, in described coupling, adopts artificial assisting sifting to combine with automatic coupling.

9. remote sensing images geometric correction method according to claim 1 is characterized in that, after described correction, corrects accuracy checking as a result, judges that whether correct the result meets the demands, if do not meet the demands, then adjusts the selection result at reference mark.

10. remote sensing images geometric correction method according to claim 1 is characterized in that, in described correction, adopts OpenMP to carry out the unit multi-core parallel concurrent at the geometric correction of piecemeal and handles.

Images