CN103218783B - Satellite remote sensing images fast geometric correcting method based on control point image database - Google Patents

Abstract

Satellite remote sensing images fast geometric correcting method based on control point image database, comprising: the geographic range of determining image to be corrected; For the geographic range of determined image to be corrected, in the image database of control point, retrieve all satisfactory control points; Carry out Auto-matching reconnaissance, be identified for the dominating pair of vertices of geometric correction; The quantity and the whether meeting geometric that distributes that judge the dominating pair of vertices of coupling are corrected requirement, if so, enter next step, if not, return to described searching step; According to sensor type and imaging characteristics, select Geometric rectification; Dominating pair of vertices based on coupling and the Geometric rectification of selection are set up the transformational relation of pixel coordinate and geodetic coordinates; Carry out geometric correction based on digital elevation model and described transformational relation, obtain the digital orthoimage through correcting. This image rectification method, is incorporated into image automatic Matching in automatically choosing of control point, can realize the automatic geometric of remote sensing image and correct processing.

Description

Satellite remote sensing images fast geometric correcting method based on control point image database

Technical field

The present invention relates to a kind of Geometric Correction of Remote Sensing Image method, more specifically relate to a kind of based on control point image databaseSatellite remote sensing images fast geometric correcting method.

Background technology

Along with the development of remote sensing technology, particularly the development of remote sensor technology, obtains by remote sensing technologyRemote sensing image or the purposes of data more and more wider. At present, the range of application of remotely-sensed data has expanded to Social Information ServiceField, for example, is widely used in mapping, agricultural, forestry, geological and mineral, the hydrology and water resource, environmental monitoring, natural calamity, districtThe aspects such as domain analysis and planning, military affairs, soil utilization. The remote sensing image with accurate geographic coding can be soil, planning, ringThe different fields such as guarantor, agricultural, forestry, ocean provide the characters of ground object and the information that need separately.

Via satellite or when the flying platform such as airborne platform obtains remote sensing image data or other data, can be subject to skyGas, daylight, the impact of external factor such as block, meanwhile, in the time of data acquisition, the height of flying platform, attitude can change,Therefore, in the time carrying out remote sensing images shooting, tend to cause the problems such as image translation, rotation, convergent-divergent. In addition, according to optical imageryPrinciple, when camera imaging, according to the imaging of central projection mode, therefore ground height rises and falls and will cause in the time of imagingThe existence of height displacement. Above-mentioned combined factors, can cause the error of remote sensing image, such as heeling error, projection error etc. Therefore,Before using these remote sensing image/data, need obtained original remote sensing image to carry out orthorectify.

Traditional remote sensing image orthorectify generally comprises: first measure by field operation or from existing topographic map dataGather ground control point (GroundControlPoint, GCP) and digital elevation model (DigitalElevationModel, DEM) information; Then these information are imported to professional remote sensing or digital Photogrammetric System; Next by loadingThe system of GCP and DEM information remote sensing image is carried out to orthorectify.

Ground control point GCP (GroundControlPoints) is that satellite remote sensing image geometry is corrected and geo-locationTime important reference data sources. Correct in processing procedure, for reaching certain correction precision, some at remotely sensing image geometricGround control point be absolutely necessary, the object coordinates by control point and corresponding picpointed coordinate build video imaging modelAnd solving model parameter or existing imaging model is optimized to compensation solves compensating parameter, finally set up object space in imaging processCorrect transformational relation with image space.

The general pattern that adopts full field measurement in control point in traditional work process, needs through " collecting existing controlData, survey that district is made an on-the-spot survey, reconnaissance is laying of markstone, field operation is measured, interior industry arrangement " etc. the operation of series of complex, although along with GPS-RTK etc.The development of advanced measuring technique, workload and complexity that field operation is measured greatly reduce, but necessary control point field measurement workDo still inevitable.

In addition, due to satellite remote-sensing image coverage large (taking No. three satellites of resource as example, coverage be 50 kilometers ×50 kilometers), want to obtain equally distributed ground control point, conventionally need to be in the scope of hundreds of kilometer and even thousands of square kilometresInterior testing, its field process amount is big, labor intensive and material resources are self-evident. In addition, at earthquake, flood, mudstoneThe uninhabited areas such as the natural calamity hotspots such as stream or virgin forest, marsh, desert, survey crew usually cannot enter realityGround is measured.

Even if control point is by manually choosing remote sensing image and topographic map same place obtains, also exist inefficiency,Precision is difficult to the problems such as guarantee. And phase or different sensors image are corrected during to areal different, there will be again heavilyThe situation of final election point.

Summary of the invention

According to embodiments of the invention, provide a kind of satellite remote sensing images based on control point image database several fastWhat correcting method. The method comprises: the geographic range of determining image to be corrected; For the geographical model of determined image to be correctedEnclose, in the image database of control point, retrieve all satisfactory control points; Carry out Auto-matching reconnaissance, be identified for how muchThe dominating pair of vertices of correcting; The quantity and the whether meeting geometric that distributes that judge the dominating pair of vertices of described coupling are corrected requirement, if so,Enter next step, if not, described searching step returned to; According to sensor type and imaging characteristics, select phaseThe Geometric rectification of answering; Dominating pair of vertices based on described coupling and the Geometric rectification of described selection set up pixel coordinate withThe transformational relation of geodetic coordinates; Transformational relation based between digital elevation model and pixel coordinate and geodetic coordinates carries out how muchCorrect, obtain the digital orthoimage through orthorectify.

According to the Geometric Correction of Remote Sensing Image method of the embodiment of the present invention, by image automatic Matching be incorporated into control point fromIn moving choosing, can realize the automatic geometric of remote sensing image and correct processing.

Brief description of the drawings

In order to be illustrated more clearly in the technical scheme of the embodiment of the present invention, below the accompanying drawing to embodiment is done to be situated between simplyContinue, apparently, the accompanying drawing in the following describes only relates to some embodiments of the present invention, but not limitation of the present invention.

Fig. 1 is according to the satellite remote sensing images fast geometric correcting method based on control point image database of the present inventionFlow chart;

Fig. 2 a-2c is exemplary control point image film and the comparison diagram of No. three satellite images of resource;

Fig. 3 is a kind of indicative flowchart of control point image film retrieval;

Fig. 4 is a kind of indicative flowchart of control point automatic matching method;

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

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

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

Fig. 8 is the schematic 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 vectogram;

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

Figure 11 is the schematic diagram that has comprised that Mismatching 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 relationship in geometric correction;

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

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

Detailed description of the invention

For making object, technical scheme and the advantage of the embodiment of the present invention clearer, below in conjunction with the embodiment of the present inventionAccompanying drawing, the technical scheme of the embodiment of the present invention is clearly and completely described. Obviously, described embodiment is thisBright a part of embodiment, instead of whole embodiment. Based on described embodiments of the invention, ordinary skillThe every other embodiment that personnel obtain under the prerequisite without creative work, belongs to the scope of protection of the invention.

Unless otherwise defined, technical term used herein or scientific terminology should be in the affiliated field of the present invention and haveThe ordinary meaning that the personage of general technical ability understands. In patent application specification of the present invention and claims, use " theOne ", " second " and similarly word do not represent any order, quantity or importance, and be just used for distinguishing differentPart. Equally, the similar words such as " " or " " do not represent restricted number yet, but represent to exist at least one.

According to the satellite remote sensing images fast geometric correcting method based on control point image database of the embodiment of the present invention,The basis that is established as with control point image database, retrieval by control point image database, control point Auto-matching withAnd target image is corrected to realize. Described control point image database is the set of control point image film. And control point imageSheet is the expansion at traditional control point, that is to say, high by orthography DOM piece and the corresponding numeral with geographical coordinate informationJourney model DEM piece replaces the data cell of single control point as database, thereby generates control point image database.

Fig. 1 schematically shows according to the flow chart of the remote sensing images geometric correction method of the embodiment of the present invention. ?Step S1, the geographic range of definite image to be corrected; At step S3, for the geographic range of determined image to be corrected,In the image database of control point, retrieve all satisfactory control points; At step S5, carry out Auto-matching reconnaissance, determine and useIn the dominating pair of vertices of geometric correction; At step S7, the quantity of the dominating pair of vertices that obtains of judgement and distribute whether meeting geometric entanglesJust requirement, if so, enters next step, if not, return to step S3 and proceed retrieval; At step S9, rootAccording to sensor type and imaging characteristics, select corresponding Geometric rectification classification; Step S11, entangles based on dominating pair of vertices, geometryPositive model is set up the transformational relation of pixel coordinate and geodetic coordinates; At step S13, based on digital elevation model (DEM) and pixelTransformational relation between coordinate and geodetic coordinates, carries out geometric correction to raw video, through the numeral of orthorectify just obtainsProjection picture.

Alternatively, before carrying out step S1, image original to be corrected is carried out to image and slightly correct. This thick correction can be adoptedBy the non-conventional method based on control point image database, for example, adopt four that in the metadata information of raw video, recordAngle point geodetic coordinates and image coordinate, build one time multinomial model, completes the geometric correction of summary, obtains not through control pointThe correction image correcting. By the scope of the thick correction raw video that effectively PREDICTIVE CONTROL point image film is corresponding, thereby dwindleThe hunting zone of Image Matching, reduces match time. Therefore, adopt the thick correction of this image can improve geometric correction of the present inventionSpeed.

According to the above-mentioned matching algorithm of the embodiment of the present invention, within the scope of certain error, can realize rapid automatized reconnaissanceOperation. Alternatively, in step S5, adopt human assistance screening to combine with Auto-matching, can improve the selection at control pointQuality, and then the precision of raising geometric correction.

In addition, alternatively, after step S13, correct result accuracy checking, judge that whether correct result meets and wantAsk, if do not met the demands, the selection result at control point is adjusted. This check result can be used as how much of correction resultPrecision evaluation reference index.

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

Control point image database builds

The fast development of modern Remote Sensing Technical, cycle and precision that remote sensing images obtain are also improved gradually, and this isThe foundation of control point image database provides reliable data base. And the development of Computer Automatic Recognition technology is controlThe effective use of point data base processed provides technical foundation. According to embodiments of the invention, set up control point image database,The attribute information at control point and image information are unified to build library management, realized and " once built storehouse, partly upgrade, repeatedly makeWith " target. , image automatic Matching is incorporated in automatically choosing of control point meanwhile, realizes the automatic of remote sensing imageOr semi-automatic geometric correction processing.

In one embodiment, each control point image film comprises two kinds of data: view data and attribute data. ItsMiddle attribute data is used for describing geographical location information, comprises four aspects:

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

2. some necessary supplementarys of geographical coordinate are described, as the coordinate system being adopted, projection pattern, ellipsoidal parameterDeng;

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

4. GCP chooses feature describes, and as the crossing of road or bridge central point, these information can be used as inquirySubsidiary conditions.

Control point image database adopts the mode of database to store, manage kimonos to a large amount of control points image filmBusiness. Traditional method is that the attribute at control point and view data are isolated and come on data storage management, in database, depositsWhat store up is the file pointer of a correspondence image, and view data is stored separately with file mode in database outside, this sideFormula has been destroyed the security of control point information integrity and database, very easily because the mistake of file is deleted and to be lost control pointImage information. Adopt the data storage shape that image and attribute are combined according to control point of the present invention image databaseFormula, adopts binary large object BLOB (BinaryLargeObject) class to image information corresponding to control point image filmType also carries out integrated storage administration as a field of list structure with attributes such as controlling coordinate, ellipsoid type, projection pattern.

To the DOM(digital orthophoto map of existing 1:10000 yardstick, DigitalOrthophotoMap), DEM(digital elevation model, DigitalElevationMode) result map arranges, and chooses that the trend of the times is strong, clean mark, spyLevy obvious area, as the areas such as intersection, bridge, ridge angle point gather shadow by a certain size (as 512 × 512 pixels)Photo obtains the information such as plane coordinates information, ellipsoid, projection from DOM, obtains respective area from corresponding dem data simultaneouslyThe height value in territory, then unifies acquisition image information and attribute information into library storage.

Fig. 2 a is that (512 × 512 pixels, resolution ratio is the control point image film cutting out from the DOM image in somewhere1 meter); Fig. 2 b be areal No. three, resource face panchromatic image (resolution ratio is 2.1 meters), the money that Fig. 2 c is arealMultispectral image (resolution ratio is 5.8 meters) is faced in source for No. three. Contrast Fig. 2 b and Fig. 2 c can choose phase very intuitively from Fig. 2 aThe same place of answering is as control point.

In addition, the coverage of control point image database and memory space are after certain scale, if physics routinelyOrder is retrieved one by one, can expend the plenty of time, is unfavorable for practical application, for quick-searching goes out control strip, in realityIn application, need to carry out partitioned storage by geographical coordinate to control point database. In subregion process, according to hypsography,The complexity of thing and concrete application, to interested area and the closeer area of distribution, control point, reduce division scope,And expand division scope for secondary sites or water field of big area etc., finally ensure that the number of control points in each region is protected substantiallyHold unanimously, accelerate the speed of retrieval.

The retrieval of control point image database

Search function is to weigh the important technology index of Database Systems. For control point image database,In the time carrying out geometric correction, can go out available control point image according to image quick-searching to be corrected be the application system constructionOne of success or failure key.

Fig. 3 is a kind of indicative flowchart of image film retrieval in control point in step S3. At step S301, based on step S1The geographical position range of middle estimated summary, carries out the retrieval of based target regional center point longitude and latitude; At step S303, according toThe attribute information of control point image film screens; At step S305, carry out content-based advanced search; At step S307,Obtain required control point image film. Introduce respectively the step of control point image film retrieval below.

First, predict the geographic range of its summary according to the orbit parameter of image to be corrected, and then estimate target areaThe general geographic location scope at interior control point. Because the general location of prediction has error, remote sensing image to be corrected exists several in additionWhat distortion, generally all can be with reference to the error radius R of an estimation in the time of the scope of the geographical position of estimation target area, and its value is generalFor 2 to 3 times of control point image size. Wherein, the upper left angle point latitude and longitude coordinates of supposing this scope is (L1, B1), the lower right cornerPoint latitude and longitude coordinates is (L2, B2), and the coordinate of target area image center to be checked is (L0, B0); L1=L0-R, L2=L0+R，B1=B0-R，B2=B0+R。

The retrieval of the center longitude of based target area image is the retrieval in the summary geographical position based on estimated.Data item L and B in the image database record sheet of control point, for the control point video search based on general location provides basicCondition. From above-mentioned coordinate relation, the expression formula for search that can set up " L1≤Li&&Li≤L2&&B2≤Bi&&Bi≤B1 " entersRow relation retrieve.

Based target image attribute retrieval, according to the resolution ratio of target image, sensor type, imaging time etc., sieveSelect available control point image film, effectively dwindle result for retrieval scope.

According to embodiments of the invention, first estimate the geographical position range of summary of remote sensing image to be corrected target area,The scope of dwindling retrieval by the retrieval of based target regional center point longitude and latitude; Then judge available control point image filmResolving range, time phase and sensor type etc., the scope of further dwindling result for retrieval. Concrete retrieval mode has: baseIn the retrieval of given coordinate range, based on sensor, time the attribute information integrated retrieval such as phase and resolution ratio, content-based (demandDistribution characteristics, color characteristic, shape facility, textural characteristics) retrieval technique. According to this retrieval flow to control point image filmAfter retrieval, the control point image quantity satisfying condition can greatly reduce, and substantially can meet the demands.

The coupling of dominating pair of vertices

In the correction procedure of remote sensing image, the step of most critical is choosing of control point of the same name, and this is to determine that image rectification certainlyThe key of movingization. The automatic correction basis of realizing image and carry out Auto-matching according to existing control point image data.

An advantage of control point image film is that existing geography information has again image texture information, therefore can adopt shadowFind the image point position of the image to be corrected that control point image film is corresponding as Auto-matching algorithm, then take out control point image filmGeographical coordinate information just can form a dominating pair of vertices that has comprised object coordinates and image space coordinate.

According to remote sensing images geometric correction method of the present invention, for control point image film and wait correct image time phase,All inconsistent features of aspect such as resolution ratio, imaging angle, have adopted how much invariant feature extractions based on SIFT and have mated calculationMethod, and adopt the methods such as Rough Fuzzy C-Mean Method to carry out automatic rejection to Mismatching point, finally adopts a classical young waiter in a wineshop or an innTake advantage of matching algorithm to carry out the coupling that becomes more meticulous of sub-pixel-level.

Idiographic flow as shown in Figure 4, at step S501, according to the metadata of the coordinate information at control point, image to be correctedInformation and imaging model etc. calculate the initial coordinate of corresponding picture point, then by the size of control point image film from shadow to be correctedOn picture, cut out image blocks to be searched; At step S503, utilize Sift algorithm treating to control point image film and after cutting to searchRope image blocks is mated, and obtains preliminary matching result information; At step S505, adopt Rough Fuzzy C-Mean Method and severalThe methods such as what constraint are carried out the rejecting of Mismatching point, retain reliably accurately matching double points; At step S507, utilize least squareAlgorithm carries out essence coupling to matching result, obtains the matching precision of sub-pixel; At step S509, by the control point that the match is successfulTo formatted output in accordance with regulations to comprising in the control point message file of controlling period, object coordinates, image space coordinate.

As shown in Figure 1, carry out above-mentioned matching operation for the multiple control points image film in image coverage to be corrected,Produce multiple dominating pair of vertices. If dominating pair of vertices quantity is inadequate, retrieves acquisition by continuation and covering image capturing range to be correctedIn more control point image film, thereby density or quantity to control point image film increases, then Auto-matchingMode obtains more dominating pair of vertices, corrects while requiring when the quantity of dominating pair of vertices and distribution meeting geometric, just can carry out oneselfThe geometric correction processing of movingization. Provide below according in the relevant concrete technology of the control point image film coupling of the embodiment of the present inventionHold.

SIFT algorithm

The thought of SIFT algorithm based on characteristics of image scale selection, sets up the multiscale space of image, under different scaleSame characteristic point detected, when determining characteristic point position, determine its place yardstick, to reach the nonshrink object of putting of yardstick,In addition, this algorithm is rejected point and the skirt response point that some contrasts are lower, and extracts invariable rotary feature descriptor to reachTo the object of anti-affine transformation. This algorithm mainly comprises: (1) sets up metric space, finds candidate point; (2) accurately determine keyPoint position, rejects point of instability; (3) determine mould and the direction of key point gradient; (4) extract feature descriptor.

1. the generation of metric space

Metric space theory its objective is the multiple dimensioned spy of simulated image data while coming across computer vision field the earliestLevy. Koendetink proves that Gaussian convolution core is unique transformation kernel of realizing change of scale in the literature, and the people such as LindebergFurther prove that Gaussian kernel is unique linear kernel.

Two-dimensional Gaussian function is defined as follows:

$G(x,y,\mathrm{\σ})=\frac{1}{{2\mathrm{\π\σ}}^2}{e}^{-(x^2+y^2)}/{2\mathrm{\σ}}^2$

σ has represented the variance of Gauss normal distribution.

One width two dimensional image, the metric space under different scale represents to 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 different scale space.

In order effectively stable key point to be detected at metric space, difference of Gaussian metric space (DoG is proposedScale-space). Utilize 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,σ)

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

Then build image pyramid, image pyramid is O group altogether, and every group has S layer, and the image of next group is by upper one group of imageDown-sampled obtaining. Fig. 6 shows one group of example of difference of Gaussian metric space (DoG) image. In Fig. 6, to metric spaceOctave, raw video, through repeatedly Gaussian convolution computing, produces the image of the metric space of a series of settings. On the rightDoG image is that the image after the gaussian filtering by closing on carries out calculus of differences generation. After every single order, Gauss's image doesThe factor be 2 down-sampled, and repeat this process.

2. spatial extrema point detects

In order to find the extreme point of metric space, whether the consecutive points comparison that each sampled point will be all with it, see itLarger or little than the consecutive points of its image area and scale domain. As shown in Figure 7, middle test point and its 8 phase with yardstick9 × 2 points totally 26 somes comparison that adjoint point and neighbouring yardstick are corresponding, to guarantee in metric space and two dimensional image space allExtreme point detected.

3. key point position is determined and point of instability rejecting

A) key point exact position is determined

Utilize Taylor's second outspread formula of metric space function D (x, y, σ) to carry out least square fitting, pass through digital simulationThe extreme value of curved surface is further determined exact position and the yardstick of key point. Coordinate and yardstick that key point is final can be as accurate asSub-pixel-level.

Launch D (x, y, σ) by Taylor's formula, sampled point initial point is:

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

To X differentiate, and to make it be zero, that is:, just can try to achieve sampling initial point position

For: $\hat{X}=-\frac{{\∂}^2{D}^{-1}}{{\∂X}^2}\frac{\∂D}{\∂X}$

Be $\left[\begin{array}{ccc}\frac{{\∂}^{2}D}{{\∂\mathrm{\σ}}^2}& \frac{{\∂}^{2}D}{\∂\mathrm{\σy}}& \frac{{\∂}^{2}D}{\∂\mathrm{\σx}}\\ \frac{{\∂}^{2}D}{\∂\mathrm{\σy}}& \frac{{\∂}^{2}D}{{\∂y}^2}& \frac{{\∂}^{2}D}{\∂\mathrm{yx}}\\ \frac{{\∂}^{2}D}{\∂\mathrm{\σx}}& \frac{{\∂}^{2}D}{\∂\mathrm{yx}}& \frac{{\∂}^{2}D}{{\∂x}^2}\end{array}\right]\left[\begin{array}{c}\mathrm{\σ}\\ y\\ x\end{array}\right]=-\left[\begin{array}{c}\frac{\∂D}{\∂\mathrm{\σ}}\\ \frac{\∂D}{\∂y}\\ \frac{\∂D}{\∂x}\end{array}\right]$

B) low contrast is rejected

ByIf, | D (X) | < 0.03, this contrast is lower, rejects.

C) removal of skirt response

An extreme value that defines bad difference of Gaussian is having larger principal curvatures across the place at edge, and is hanging downThe direction of straight edge has less principal curvatures. 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 obtained by the adjacent poor estimation of sampled point.

The principal curvatures of D and the characteristic value of H are directly proportional, and making α is eigenvalue of maximum, and β is minimum characteristic value:

Tr(H)=D_xx+D_yy=α+β

Det(H)=D_xxD_yy-(D_xy)²=αβ

Make α=γ β:

$\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 characteristic values are equal, increase along with the increase of r. Therefore, in order to examineSurvey 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

The gradient direction distribution characteristic of utilizing key point neighborhood territory pixel is each key point assigned direction parameter, makes operator toolStandby 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$

Above formula is mould value and the direction formula that (x, y) locates gradient. Wherein L yardstick used is each key point place separatelyYardstick.

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 for working asThe position of front key point. First, reference axis is rotated to be to the direction of key point, to guarantee rotational invariance. Next with keyCentered by point, get 8 × 8 window. In Fig. 8, each little lattice represent a pixel of key point neighborhood place metric space, arrowDirection represents the gradient direction of this pixel, and arrow length represents gradient-norm value, and the scope that in figure, circle represents Gauss's weighting (morePixel gradient directional information contribution near key point is larger). Then on every 4 × 4 fritter, calculate the gradient side of 8 directionsTo histogram, draw the accumulated value of each gradient direction, can form a Seed Points, as shown in Fig. 8 right half. In this figure oneIndividual key point by 2 × 2 totally 4 Seed Points form, each Seed Points has 8 direction vector information. This neighborhood directivity informationThe thought of associating has strengthened the antimierophonic ability of algorithm, also provides good for the characteristic matching that contains position error simultaneouslyFault-tolerance.

In computational process, in order to strengthen the robustness of coupling, alternatively, to each key point use 4 × 4 totally 16 kindsSon point is described, and just can produce 128 data like this for a key point, the SIFT feature that finally forms 128 dimensions toAmount. Now SIFT characteristic vector has been removed the impact of the geometry deformation such as dimensional variation, rotation factor, then continue by feature toThe length normalization method of amount, the impact that can further remove illumination variation.

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

Based on the Remote Sensing Images Matching of SIFT descriptor

After the piece image to independent carries out above-mentioned feature extraction and feature description, just obtain in this figure allFeature and descriptor thereof, being made as image 1(is realtime graphic), its characteristic point quantity is m. Accomplish two width images match to riseCome, namely obtain the pixel (referring to characteristic point here) mating in two width images. First, be to other piece image (image2, i.e. reference picture) carry out identical feature extraction and feature is described process, obtain the feature that quantity is n; Secondly be exactlyWill be in image 1 in n characteristic point, search out the characteristic point of r to coupling in m characteristic point and image 2, wherein r≤m, r≤n,And in order to calculate the geometrical relationship between image 1 and image 2, should ensure r >=8 robust, if can not findTo the r satisfying condition, just need to adjust the precision of coupling, such as the similitude requirement reducing between matching characteristic point, make energyObtain more matching characteristic point right.

The coupling of characteristic point is in fact the coupling of its descriptor, and the descriptor of characteristic point is that this feature has been carried out in factA quantitative description, can be applied to matching algorithm. The coupling of characteristic point descriptor is to carry out in description space in fact, such as SIFT Feature Descriptor is the vector of 128 dimensions, therefore SIFT descriptor mates in 128 dimension spaces. ?Describe in space, the matching degree of characteristic point descriptor is measured with distance, and two nearest descriptors are generalJust represent a pair of characteristic point of coupling. And in description space, generally have the definition of following two kinds of distances:

A) Euclidean (Euclidean) distance. Be p dimensional vector for 2 x in p dimension space and y(), there are theirsEuclidean distance is 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 formula₁,…,x_p)^T,y=(y₁,…,y_p)^T. Euclidean distance can be applicable to, based on histogrammic descriptor, retouchEvery one dimension of stating son has identical weight, such as SIFT descriptor, and GLOH descriptor and PCA-SIFT descriptor etc.

B) Ma Shi (Mahalanobis) distance. If every one dimension of descriptor has different weights, need to use Ma ShiThe distance of range measurement between them, the weight vectors of now establishing descriptor is s=(s₁,…,s_p)^T, 2 points in P dimension spaceMahalanobis distance between x and y is:

$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 formulaMahalanobis distance can be applied to the controlled filtering descriptor of direction, the constant description of differenceSon, the constant descriptor of square and complex filter descriptor etc.

In an embodiment according to the present invention, use SIFT Feature Descriptor to mate, it is using Euclidean distance as commentingSentence the standard of descriptor matching degree. Figure 10 shows the coupling based on Sift descriptor according to an embodiment of the invention and calculatesMethod flow process. Part in Figure 10 in square frame represents that the extraction of SIFT characteristic point and feature describe process.

According to the algorithm shown in Figure 10, before output matching result also will by " minimum distance with second closely thanCompared with back-and-forth method " select match point. Particularly, first set a rational threshold value t, then relatively in A and n point, ownThe distance of point, then finds out with A apart from minimum some B and the some C little with A distance second, and if d_ABd_ACWhen < t,Just think that A and B are reasonable correct match points. The benefit of doing is like this in the time that A has a lot of similar coupling, n put in AMinimum distance and differing while not being a lot of with the little distance of A second, can think that this is not a reasonably coupling. Only having ought beLittle distance time mutually far short of what is expected, when the minimum point of distance " holds a safe lead " other, is just admitted it with the second little distance, because ofFor this basic guarantee this be a highly stable correct coupling instead of an ambiguous coupling.

Mismatching point is rejected

To rear, substantially realize images match at the matching characteristic point that has obtained realtime graphic and reference pictureTask, but the characteristic point of coupling can only representative image between local relation, the matching characteristic point of limited quantity can notReflect completely image between holotopy. In corresponding matching double points, may there is mistake coupling or matching errorLarge point is right, due to these put on existence will affect matching precision, for the application based on match point file impacts, soAfter initial matching finishes, explore the right elimination method of Mismatching point also extremely important. Further, Mismatching point automatic rejection energyEffectively automation ground improves the accuracy rate of coupling.

Figure 11 shows the matching algorithm flow process that has comprised according to an embodiment of the invention Mismatching point rejecting step.In Figure 11, first carry out Feature Points Matching according to aforementioned SIFT algorithm based on how much invariant features, then for matching resultScreen deletion, wherein relate to Rough 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 SIFT algorithm to extract respectively real-time imaging and the characteristic vector information with reference to image. According to twoThe SIFT characteristic vector of width image, adopts the Euclidean distance of key point characteristic vector to be used as the similar of key point in two width imagesSex determination tolerance, Specific Principles is that arest neighbors (NN) and the second neighbour's (SCN) ratio of distances constant (NN/SCN) is minimum as two widthThe similarity determination of key point tolerance in image. Afterwards, utilize Rough Fuzzy C-Mean Method to carry out tentatively picking of Mismatching pointRemove.

The Rough Fuzzy C-Mean Method that Mismatching point is rejected

The basic thought of cluster analysis is very simple, directly perceived and simple, and it is the pattern feature phase to be sorted according to eachClassify like degree, the similar class that is classified as, dissimilar as an other class. Cluster analysis comprise two substantially inHold: the tolerance of pattern similarity and clustering algorithm. FCM (FuzzyC-Means) algorithm is a kind of with flat in infima speciesSide's error and be clustering criteria, calculates the degree of membership that each sample belongs to each fuzzy subset (cluster), minimum by object functionPickard iteration between the necessary condition of changing is carried out implementation algorithm.

Definition 1: establish X={x₁,x₂,x₃,…，x_nBe object set to be sorted, for i class w_i, its barycenter is v_i, definition For can certainly belonging to w_iThe object set of class:

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

B) simultaneously in order to ensure A_iBe unlikely to obtain too small, need meet for, at least exist k ∈ 1,2 ..., n}, makes

Wherein A_iBecome upper approximate limit, upper approximate limit has been portrayed the border that likely belongs to the object of i class, if certainObject does not belong to the scope that approximate limit defines, and it belongs to the negative territory of this class, does not belong to this class completely.

Definition 2: the object function of Rough Fuzzy C-average (RoughFuzzyC-Means, RFCM) 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{ij}}^2$

Constraints 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 above formula is:

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

Above formula be can be calculated:

$u=1/\underset{k=1}{\overset{c}{\mathrm{\&Sigma;}}}{(d_{\mathrm{ij}}^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 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 the object that belongs to certain class has been divided into sure, possible and three of negatingSet, carries out cluster taking square error in the infima species of all possible object with as clustering criteria. RFCM algorithm and FCM algorithmMaximum difference is, it thinks x_jBelong to w_iDegree of membership u_ijCalculating only comprise x with upper in approximate_jClass relevant, if certainClass w_iUpper approximate in do not comprise x_j, this class is to x_jDegree of membership without any contribution.

The object function of RFCM algorithm be find in each cluster infima species square distance and, it has been got rid of beyond suprasphereThe object that can not belong to class. For u_ijCalculating, if for class w₁, object x_jOnly be positioned at the upper approximate of class, in,?, i.e. u_1j=1, if x_jBelong toWithCommon factor, u_1j，u_2jCalculating only with class w₁, and w₂Relevant, andWith w₃Irrelevant, that is:

$u_{1j}=1/{(1+\frac{{\left|\right|x_j-v_1\left|\right|}^2}{{\left|\right|x_j-v_2\left|\right|}^2})}^{\frac{1}{m-1}},$ $u_{1j}=1/{(1+\frac{{\left|\right|x_j-v_2\left|\right|}^2}{{\left|\right|x_j-v_1\left|\right|}^2})}^{\frac{1}{m-1}}$

If x_jDo not existIn, u_1j=0。

Object function, provided u_ijMore new formula. According to this formula, with v_jCentered by u in the radius suprasphere that is α_ijBe 1.

Based on above-mentioned mathematical analysis, according to the Rough Fuzzy C-Mean Method mistake under the geometrical constraint of the embodiment of the present inventionJoining a rejecting step mainly comprises:

According to the matching double points generating and the resolution information of image, calculate the geometrical constraint condition between image:

● in resolution ratio situation, ask slope k that corresponding points are right and apart from S;

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

Then, utilize Rough Fuzzy C-Mean Method (RFCM), ask the degree of membership of (k, S) or (X, Y), to all couplingsPoint, to carrying out cluster analysis, is deleted inhomogeneous point right, only retains the point being included in sure gathering set.

I. definite class is counted c (2≤c≤N), parameter m, initial matrix, the upper approximate boundaries A of class_iWith a suitable decimal, s=0；

Ii. calculate barycenter

If iii.U_ij=0, otherwise upgrade

If iv. || U^s-U^s+1|| < ε, stop, otherwise s=s+1, forwards b to.

As shown in figure 11, if still have and exceed setting number after Rough Fuzzy C-Mean Method is carried out match point examinationThe match point of (for example, 7), adopts random sampling consistency detecting method to carry out further match point and screens out.

Random sampling consistency detecting method (RANSAC method)

A) amount of calculation of RANSAC algorithm

In RANSAC algorithm, require to ensure under certain fiducial probability in the sampling of M group, to have the data of one group of sampling at leastInterior point (inliers) entirely. Utilize the following formula can be in the hope of the minimum sampling number M meeting the demands.

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

Wherein: ε is data error rate (exterior point (outliers) is in the shared ratio of initial data), m is computation model ginsengThe minimum data amount that number needs, P is fiducial probability. As can be seen from the above equation, M and ε, m, P is exponential relationship. Following table shows P=0.95 o'clock, M was with m, the situation that ε changes.

As can be seen from the above table, in the time that model more complicated, ε are higher, M is very large, directly causes under RANSAC efficiency of algorithmFall. Need the time if randomly draw one group of sampling from initial data; Need time t by one group of sample calculation model parameter; WithThe data testing model parameter average out to t that takes time, uses N data detection (all data inspection) to need time N_t. ThereforeRANSAC algorithm is required computing time:

T=M(t_s+t_c)+MN_t

Wherein, M (t_s+t_c) be the time that the sampling of M group is extracted and model parameter calculating needs, MN_tFor M model parameter inspectionTest the time needing.

Find out from this formula, the time that 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 model parameter inspection needs. The in the situation that, data error rate definite at model being definite,

In order to ensure the fiducial probability of result, M can not reduce.

Therefore,, in order to improve the efficiency of algorithm, can only participate in from reducing model parameter quantity, the minimizing model parameter of inspectionThe time that inspection needs sets out.

B) RANSAC algorithm steps

Calculate minimum sampling number M according to fiducial probability P and data error rate ε; Calculate model parameter corresponding to sampling, use instituteThere is initial data testing model parameter quality, obtain the inliers quantity of each model parameter; According to inliers quantity and mistakePoor variance is selected optimal models parameter; The inliers corresponding by optimal models parameter estimates final mask parameter.

Least square method essence coupling

Carry out further after match point selects through RANSAC algorithm, can also utilize young waiter in a wineshop or an inn for fitting of a polynomialTake advantage of method deletion fitting residual error to be greater thanThe matching double points of error doubly.

As shown in figure 11, less if coupling is counted, be for example not more than 7, alternatively, can adopt parameter adaptive SIFTAlgorithm increases the match point for screening.

Satellite remote-sensing image is because extraneous factors such as being subject to weather, sunlight, block has a strong impact on, and existence is because of differenceThe factor such as imaging time, angle, distance and between the problem such as the image translation that causes, rotation, convergent-divergent and various sensorThe difference of imaging pattern camera parameter model trajectory, even if adopt SIFT algorithm, under extreme case, also may occur extractingThe few or basic extraction of characteristic point less than the phenomenon of characteristic point, cause mistake coupling or it fails to match. For this situation, can adoptWith improved parameter adaptive SIFT image matching method, according to different imaging characteristics and quality, choose corresponding strategy alsoDetermine the respective threshold of feature point detection according to relevant weight definition rule.

The parameter relating in SIFT feature extraction and matching process is more, for verify each parameter feature extraction andJoin the effect of playing in process, carried out Matching Experiment, main by the size of certain parameter to be verified of amendment, to this ginsengNumber regulates according to theoretical rule change in algorithm, meanwhile keeps other parameter constants, each by coupling counting statisticsThe match point number in parameter when size, Mismatching point number, match time etc., and carry out statistical analysis.

By known to the analysis of statistical results of experiment, its universal law is as follows:

● in the time that the sampling interval of every rank (Octave) increases gradually, match point number increases rapidly, Mismatching point quantityDo not change, increase match time;

● increase gradually Gaussian convolution core σ, mate the entirety of counting and become downward trend, but result is best in the time of σ=1.6;

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

● along with the increase gradually of bent ratio threshold value, match point number increases gradually, and Mismatching point number is almost constant;

● coupling threshold value minimum distance increases with the ratio of time minimum distance, and match point number increases gradually, mistakeJoining counts also progressively increases.

Based on above-mentioned experimental result, parameter adaptive SIFT algorithm is proposed. The step of parameter adaptive SIFT algorithm isInformation by image self is calculated the threshold value of automatically determining each parameter, thereby plays by changing parameter the feature that increase is extractedThe object of some quantity, increase match point quantity. Specifically can comprise the steps:

A) calculate input image average gray or image is done to Auto Laves, if then by average gray divided by 10Result is less than 1.0, reduces the threshold value of contrast, increases the characteristic point quantity of extracting;

B) in the situation that meeting matching precision, if counting, coupling is less than at 7 o'clock, increase bent ratio threshold value;

If c) through a) b) feature quantity that two steps are extracted when still less or number of matches is less, do not affecting coupling essenceIn the situation of degree, suitably amplify the proportion threshold value of minimum distance and time minimum distance.

As shown in Figure 1, carrying out, after the coupling of control point, according to sensor type and imaging characteristics, selecting several accordinglyWhat corrects model (step 9), and then the dominating pair of vertices based on mated and selected Geometric rectification are set up pixel coordinateAnd the transformational relation (step 11) between geodetic coordinates, finally based on digital elevation model (DEM) and pixel coordinate and the earth seatTransformational relation between mark carries out geometric correction (step 13) to raw video. The selection of Geometric rectification is mainly described belowStep with geometric correction.

Geometric rectification and selection thereof

Linear array push-broom type imaging sensor be line by line with time sequential mode obtain two dimensional image. It is generally elder generation's shape in image planesImage into a line, then satellite is pushed ahead along predetermined track, forms one by one a width bidimensional image, imaging side after scanningFormula as shown in figure 12. On image every a line pixel synchronization imaging and centered by projection, whole image is multicenter projection.

In Figure 12: p_kFor arbitrary picture point on image, x_kFor the x coordinate of the upper imaging point of scan line k, c is the main distance of sensor,O_kFor the projection centre of scan line k, o_kFor the principal point of scan line k, l_kFor scan line k is from projection centre O_kThe light sending.

In geometric correction, a problem the most basic will be set up rational remote sensing image imaging model exactly, so-called remote sensing shadowThe imaging model of picture refers between the topocentric geodetic coordinates (X, Y, Z) that the coordinate (x, y) set up on image is corresponding with itMathematical relationship. That is:

X=f_x(x,y,g)

Y=f_y(x,y,h)

Wherein g, the impact that h is other factors.

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

Physical model refer to be considered to as time cause the physical significance of deformation of image as surface relief, earth curvature, atmosphereThe factors such as refraction, camera distortion, then utilize these physical conditions to be construed as picture geometrical model, and the most representational is conllinear barPart equation is the strict imaging model on basis.

Universal model is not considered imaging mechanism, but directly describes the geometrical relationship of picture point and object point with mathematical function,It has the feature of generality, confidentiality, high efficiency. Conventional universal imaging model has multinomial, direct linear transformation, affineConversion, rational function model etc.;

In an embodiment according to the present invention, taking strict imaging model, rational function model and multinomial model as how muchThe optional model of correcting. According to sensor type and imaging characteristics, select corresponding Geometric rectification, for example, SPOT5 satelliteImage, generally selects strict imaging model, selects rational function mould for satellite images such as IKONOS, P5, RapidEyeType, nothing money rail parameter or rational function model coefficient file, multinomial model selected.

Picture point and object space point strict geometric process one to one when what strict imaging model was described is imaging, with what commonly useCollinearity condition equation is representative. According to one embodiment of present invention, the high-precision Satellite providing according to SPOT5 is counted one by oneAccording to, attitude data and time data, by imaging process, set up picture point (l, p) and ground point (X, Y, Z) orBetweenMathematical relationship. Figure 13 shows picture point and topocentric mathematical relationship.

Strict imaging model coordinate system used 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 to be determined by the attitude control system (ACS) of satellite,Wherein initial point is at centroid of satellite, X₁Axle, Y₁Axle, Z₁Axle is got respectively three principal moments axles of satellite. X₁Axle is along satellite transverse axis, Y₁AxlePoint to satellite flight direction (not exclusively overlapping), Z along the longitudinal axis₁Axle is determined according to right-handed system rule. Satellite attitude measurement is at thisIn body coordinate system, carry out, three parameters describing its spatial attitude are: pitching (pitch) is around body coordinate system X₁Revolving of axleTurn, 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)_pWhile representing respectively P components and parts imaging of charge-coupled device (CCD), direction along ng a path is with verticalThe inclination angle of orbital direction.

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

O₂-X₂Y₂Z₂The direction of each axle is the position in the t moment by satelliteAnd speedTwo direction vectors are trueFixed, 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 are to be all defined in earth agreement reference frame (ITRS),ITRF is the realization of ITRS, is to be responsible for setting up and safeguarding by the Terrestrial Reference Frame portion of IERS central office, this reference frameStation coordinates and speed are to utilize the space geodesy technique such as to defend such as very long baseline interferometry(VLBI (VLBI), Laser Measuring and gatherData calculate, be a reference frame constantly changing, that adopt at present is ITRF90, calculates one for convenientAs adopt WGS84 coordinate system to replace the definite coordinate system of ITRF framework. This be because, WGS84 is the earth's core body-fixed coordinate system,ITRF is inertial coodinate system, and what final mapping needed is body-fixed coordinate system, and ITRS and WGS84 only differ 1 meter of left and right, fromITRF90 is as shown in the table to the conversion parameter of WGS84:

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\\ T_3\end{array}\right)+(1+D)\left(\begin{array}{ccc}0& {-R}_3& R_2\\ R_3& 0& -R_1\\ -R_2& R_1& 0\end{array}\right)\left(\begin{array}{c}{X}_{\mathrm{ITRF}}\\ Y_{\mathrm{ITRF}}\\ Z_{\mathrm{ITRF}}\end{array}\right)$

WGS84 coordinate system be one by global the earth's core coordinate frame of reference and one group of corresponding model and WGS84 the earthThe measurement reference system that level surface forms. The initial point of this coordinate system is positioned at earth centroid, and Z axis is identical with reference to utmost point sensing with IERS, XAxle points to IERS guide meridian and crosses initial point and the intersection point perpendicular to Z axis plane, and Y-axis is determined by right-handed system rule;

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

I, the corresponding imaging time of calculating l capable p row pixel

t=t_c+lsp×(l-l_c)

Wherein, t_cFor the photography time of central scan row, the sweep time that lsp is every row, l_cFor the line number of central scan row

II, interpolation go out the position of t moment satellite

$\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)}$

III, the imaging inclination angle (ψ being listed as by p_x)_p(ψ_y)_pCalculate (l, p) 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 IV, 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,WithRespectively that the position of pixel (l, p) imaging moment t satellite and velocity interpolation are while going out tCarve the attitude angle of satellite:

$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 moment satellite can be obtained from auxiliary data file.

By the imaging direction vector under navigation coordinate systemBe converted into the imaging direction vector under 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 V, 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₂For orbital coordinate system O-X₂Y₂Z₂Three change in coordinate axis direction vector VI, topocentric location:

As shown in figure 15, measurement vector u₃With the intersection point on ground be M, the position vector of satellite isM is connected with the earth's coreCorresponding vector is. Again because M point is a+h at major axis, on the ellipsoid that minor axis is b+h, so have after simultaneous equations:

$\stackrel{\&RightArrow;}{O_{3}M}=\stackrel{\&RightArrow;}{P}\left(t\right)+\mathrm{\&mu;}\&times;u_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{X^2+Y^2}{A^2}+\frac{Z^2}{B^2}=1$

Wherein,

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

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

$u_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 acquisition navigation coordinate system and the observed direction of earth axes; Projection centre is(X_P,Y_P,Z_P), coordinate vector is (X, Y, Z), and projection coefficient is u, and spin 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 obtain 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 thus the topocentric the earth corresponding with it of the coordinate (x, y) on image in remote sensing image imaging processMathematics mapping relations between coordinate (X, Y, Z).

RFM model (RationalFunctionModel, rational function model) is one in general purpose transducer modelKind, its essence is to describe strict imaging process by more complicated Mathematical Modeling, sets up object space coordinate and image space pictureThe direct corresponding relation of point coordinates. RFM model can effectively seal the concerning security matters such as satellite orbit parameter, attribute sensor parameterDress, user can be under the prerequisite of not knowing relevant information, directly just can calculate the space of picture point according to RFM model and parameterCoordinate, thus carry out various application. RFM model generally adopts RPC parameter (RationalPolynomialCoefficients, rational polynominal coefficient) represent, RPC parameter generally has 80.

A, the RFM modular concept of utilizing RPC to locate

RFM model is picpointed coordinate (s, l) to be expressed as to (wherein φ is for warp with millet cake space coordinates (φ, λ, h) correspondinglyDegree, λ is latitude, h is elevation) be the polynomial ratio of independent variable:

$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) is the ground point space coordinates (φ, λ, h) after standardization, and (s, l) is the picpointed coordinate after standardization(S,L)，a_i,b_i,c_i,d_i(i=1～20),φ₀,φ_s,λ₀,λ_s,h₀,h_s,S₀,S_s,L₀,L_sFor the RPC parameter of photograph, φ₀,λ₀,h₀，S₀，L₀For standardization translation parameters, φ_s,λ_s,h_s,S_sL_sFor standardization scale parameter. They are carried together with image dataSupply.

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), Denl (U, V, W) can be expressed as three times about (U, V, W) tooPolynomial form.

Solving of b, RPC parameter

By 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

Error equation is:

V=Bx-l,W

In above 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 weight matrix.

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

(B^TB)x=B^Tl

When the employing second order for resolving RPC model parameter or second order are when above, there is overparameterization in RPC model,In RPC model, the variation of denominator is very violent, so just causes design matrix (B^TB) state variation, design matrix becomes veryDifferent matrix, causes least square adjustment not restrain.

The method that adopts Spectrum correction iteration solving method equation, no matter normal equation is good state, morbid state or rank defect, it resolves journeyOrder does not all need to add any in the time that normal equation is good state, just can converge to exact solution through iteration several times. In the time of normal equation wan, receiveThe speed of holding back is slightly slow, and estimated result is without partially. Concrete formulae express is:

Be provided with (B^TB)x=B^TL, formula both sides add x simultaneously above, obtain

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

Wherein, E is t rank unit matrixs, because unknown parameter x is all contained on formula both sides, so can only adopt the side of iterationMethod solves, and 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, directly adopts a suitable mathematics multinomial of number of timesCoordinate relation between image respective point before and after correcting is described.

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 formula, (u is v) the picture plane coordinates of picture point; (X, Y) is its geodetic coordinates of millet cake accordingly; a_ij,b_ijFor manyThe coefficient of item formula. Here polynomial exponent number is generally not more than three times, because more the multinomial of high-order often can not improve precisionCan cause on the contrary the relevant of parameter, cause the reduction of model orientation precision.

The conversion of setting up pixel coordinate and geodetic coordinates at the Geometric rectification of the dominating pair of vertices based on coupling and selectionAfter relation, at step S13, the transformational relation based between digital elevation model and pixel coordinate and geodetic coordinates carries out geometry and entanglesJust, obtain the digital orthoimage through orthorectify. The elevation model is here to obtain from existing basic surveying and mapping product figureArrive.

In addition, alternatively, after step S13, correct result accuracy checking, judge that whether correct result meets and wantAsk, if do not met the demands, the distribution to control point and quantity selection result are carried out manual intervention adjustment.

According to embodiments of the invention, state in realization satellite remote sensing images based on control point image database several fastWhen what correcting method, can adopt the parallel accelerated processing method of geometric correction based on OpenMP.

The parallel acceleration of geometric correction based on OpenMP processed

The specification of OpenMP is initiated by SGI, and it is a kind of multiprocessing towards shared drive and distributed shared memoryDevice multithreaded programming language. OpenMP has good portability, supports Fortran and C/C++ programming language, operationUnix system and Windows system are supported in system platform aspect. The importance of OpenMP is, it can be multi-thread for writingCheng Chengxu provides a kind of simple method, carries out complicated thread creation, synchronous, load balance and destruction work without programmerDo. OpenMP is for For Do statement particularly suitable, and it can cross by add several simple statements in original program just realShow the parallel processing of multinuclear, make full use of the cpu resource of computer, realize the acceleration of calculating. Therefore, for the geometric correction of piecemealIt will be very effective adopting OpenMP to carry out the processing of unit multi-core parallel concurrent, and hardware is not had to too much requirement yet.

Remote sensing image is actually from a geometric space and is transformed into another geometric space at correction procedure, once set upCorrect model, this conversion is unique, only relevant with position, and in addition because imaging is continuous, after correcting, image is alsoContinuous, therefore can, by the thought of piecemeal, by calculating four coordinates before and after angle point conversion, then set up in piece relativeSimple affine Transform Model, thus the calculating of correcting model simplified, reduce amount of calculation. Consider the impact of landform, piecemealSize is to need to consider, should be that the local piecemeal that hypsography is large is little in theory, and the little local piecemeal that rises and falls is large, and such asThe resolution ratio such as SPOT5, P5 is all 2.5 meters of left and right, thus finally divide block size to be decided to be 15 × 15, and large for HJ_1A 1B piecemealLittlely be decided to be 5 × 5. On this basis, in order to consider the needs to large image processing, can carry out piecemeal behaviour to original imageDo, be divided into some by original image, each piece is regarded piece image as and is corrected separately and then write disk, finally synthetic.

Following table shows the example that carries out the geometric correction processing based on OpenMP for different images type.

As can be seen from the above table, for strict imaging model, adopt piecemeal to add parallel processing and can significantly improve calculating effectRate, speed-up ratio is in 50 left and right, and for rational function model and its speed-up ratio of multinomial also in 4 left and right.

According to the Geometric Correction of Remote Sensing Image method of the embodiment of the present invention, by image automatic Matching be incorporated into control point fromIn moving choosing, can realize the automatic or semi-automatic geometric correction processing of remote sensing image. Particularly, the method is based on control pointImage database, obtains control point coupling number by the query and search for image database for control point, control point image Auto-matchingAccording to; Recycle the geometric correction of the multiple models such as strict imaging model, rational function model, multinomial model; Realize completeRemote sensing image automatic geometric correct solution, effectively utilized existing control point resource, improved production efficiency. WithTime, the method can be calculated fast by unit multi-core parallel concurrent, thereby further improves the correction efficiency of remote sensing image.

The above is only exemplary embodiment of the present invention, but not for limiting the scope of the invention, thisBright protection domain is determined by appended claim.

Claims (8)

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

Determine the geographic range of image to be corrected;

For the geographic range of determined image to be corrected, in the image database of control point, retrieve all satisfactory controlsPoint processed;

Carry out Auto-matching reconnaissance, be identified for the dominating pair of vertices of geometric correction;

The quantity and the whether meeting geometric that distributes that judge the dominating pair of vertices of described coupling are corrected requirement, if so, enter nextIndividual step, if not, described searching step returned to;

According to sensor type and imaging characteristics, select corresponding Geometric rectification;

Dominating pair of vertices based on described coupling and the Geometric rectification of described selection are set up turning of pixel coordinate and geodetic coordinatesChange relation;

Transformational relation based between digital elevation model and described pixel coordinate and geodetic coordinates carries out geometric correction, obtain throughCross the digital orthoimage of orthorectify,

Wherein, select strict imaging model for SPOT5 satellite image, the high-precision Satellite providing according to SPOT5 is counted one by oneAccording to, attitude data and time data, by imaging process, set up picpointed coordinate and WGS84 coordinate under ITRF90 coordinate systemMathematics transformational relation between topocentric coordinates under system, based on described digital elevation model and described picpointed coordinate and described groundDescribed mathematics transformational relation between millet cake coordinate, comes described by the unit multi-core parallel concurrent piecemeal processing based on OpenMPSPOT5 satellite image carries out geometric correction,

Wherein saidly by the unit multi-core parallel concurrent piecemeal processing based on OpenMP, described SPOT5 satellite image is carried out how muchCorrection comprises: taking 15 pixel × 15 pixels as unit, described SPOT5 satellite image is carried out to piecemeal, by calculating each image blockCoordinate before and after four angle point conversions, sets up the affine Transform Model for this image block, concurrently to image block described in eachCorrect separately, synthetic more afterwards.

2. remote sensing images fast geometric correcting method according to claim 1, is characterized in that, described retrieval comprises:

Estimate the geographical position range of summary of image to be corrected;

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

Screen according to the attribute information of control point image film;

Carry out content-based advanced search.

3. remote sensing images fast geometric correcting method according to claim 2, is characterized in that, described according to control point shadowThe attribute information of photo screen comprise according to the resolution ratio of target image, sensor type, imaging time screen can profitWith control point image film.

4. remote sensing images fast geometric correcting method according to claim 2, is characterized in that, described content-based heightLevel retrieval comprises based on demand distribution characteristics, color characteristic, shape facility, textural characteristics to be examined control point image filmRope.

5. remote sensing images fast geometric correcting method according to claim 1, is characterized in that, described Auto-matching bagDraw together:

According to the metadata information of the coordinate information at control point, image to be corrected and imaging model calculate corresponding picture point at the beginning ofBeginning coordinate, then cuts out image blocks to be searched by the size of control point image film from image to be corrected;

Image blocks to be searched after utilizing Sift algorithm to control point image film and cutting is mated, and obtains preliminary couplingObject information;

Adopt Rough Fuzzy C-Mean Method and geometrical constraint method to carry out the rejecting of Mismatching point, retain reliably accurately couplingPoint is right;

Utilize least-squares algorithm to carry out essence coupling to matching result, obtain the matching precision of sub-pixel;

By the dominating pair of vertices that the match is successful formatted output in accordance with regulations to comprising the control of controlling period, object coordinates, image space coordinateIn dot information file processed.

6. remote sensing images fast geometric correcting method according to claim 5, is characterized in that, the described Sift algorithm that utilizesMate and comprise:

Set up metric space, find candidate point;

Accurately determine key point position, reject point of instability;

Determine mould and the direction of key point gradient;

Extract feature descriptor.

7. remote sensing images fast geometric correcting method according to claim 1, is characterized in that, in described coupling, adoptsCombine with Auto-matching with artificial assisting sifting.

8. remote sensing images fast geometric correcting method according to claim 1, is characterized in that, after described correction,Correct result accuracy checking, judge whether correct result meets the demands, if do not met the demands, the selection to control pointResult is adjusted.