CN105160684A - Online automatic matching method for geometric correction of remote sensing image - Google Patents

Abstract

An online automatic matching method for geometric correction of a remote sensing image relates to the technical field of remote sensing image matching. The method has the advantage that uniformly distributed control points are efficiently and reliably obtained by virtue of free or low-cost network image map resources by fully utilizing prior geometric information of the remote sensing image. The method comprises: dividing a to-be-matched image into a plurality of regions according to the number of control points required to be acquired by utilizing an initial imaging model of the remote sensing image, and taking an image block with a certain size in each region as a processing unit; determining an approximate range of a reference image block according to the range of a to-be-matched image unit and the initial imaging model, then downloading the reference image block through a network and re-sampling the reference image block to ensure that the resolution ratio of the reference image block is close to that of the to-be-matched image unit; and performing SIFT matching on the to-be-matched image unit and the reference image block, removing gross error points, finally selecting an optimal matching point from the rest of matching points, performing least square matching, and correcting point coordinates.

Description

A kind of on-line automatic matching process for remotely sensing image geometric correction

Technical field

The present invention relates to a kind of remote sensing image automatic matching method of practicality, network Aeronautics and Astronautics photomap can be utilized to carry out on-line automatic coupling to remote sensing image.Can be applicable to the fields such as remote sensing, photogrammetric, mapping, image procossing.

Background technology

Image Auto-matching and geometry correction are the key links in Photogrammetry and Remote Sensing task, and they are visual fusion, inlay, change the basis of the senior application such as detection, map rejuvenation.Although very many to the research of image Auto-matching in the past few decades, the Auto-matching of remote sensing image is still very challenging.The performance that practical automatic matching method all should have in efficiency, robustness, precision etc., but because remote sensing image data amount is large, scene large, obtain the features such as condition is changeable, geometry deformation is complicated, existing method is difficult to take into account this several respects.In addition, the preparation with reference to image is also a difficult point of remote sensing image Auto-matching and geometry correction, and especially high resolving power often needs very high cost with reference to the acquisition of image.

In view of the limitation of existing automatic matching method, a kind of on-line automatic matching process for remotely sensing image geometric correction has important practical value.

Summary of the invention

The object of the invention is to solve the deficiencies in the prior art, a kind of quick, sane, accurate on-line automatic matching process for remotely sensing image geometric correction is proposed, the method can process the remote sensing image of arbitrary size, utilize network video map to be automatically matched to some, well-distributed high-precision control points in the short period of time, be directly used in the geometric accurate correction of remote sensing image.The advantage of the method is mainly the priori geological information making full use of remote sensing image, and by network video map resource that is free or low cost, high efficient and reliable ground obtains equally distributed reference mark.

For solving the problem, the invention provides a kind of on-line automatic matching process for remotely sensing image geometric correction, the method comprising the steps of:

Remote sensing image to be matched is evenly divided into several regions by the number at the reference mark S1. gathered as required;

If S2. all imagery zones are all processed, then complete whole Image Matching process; Otherwise, start to process next imagery zone;

S3. current image region is divided into several image units by the size of 256 pixel × 256 pixels;

If all image units S4. in current image region are all processed, then mark this imagery zone and complete process, and proceed to step S2; Otherwise, start to process next image unit;

S5: according to the approximate range with reference to image blocks in the scope of current image unit to be matched and initial imaging model computational grid photomap, then by reference image blocks corresponding to web download and be the resolution close with image unit to be matched by its resampling;

S6: utilize SIFT to mate operator and mate to image unit to be matched with reference to image blocks, and excluding gross error point, if this step obtains the match point of more than 4, then carry out step S7, otherwise proceed to step S4;

S7. from the match point that step S6 obtains, select optimal matching points and Least squares matching is carried out to it, revising the position coordinate of SIFT feature point, this is added in matching result match point.

Wherein, step S5 comprises further:

S5.1 closes on level of zoom most according to the resolution computational grid photomap of image unit to be matched;

S5.2 calculates the corresponding width with reference to image blocks and height;

S5.3 calculates the latitude and longitude coordinates of the corresponding central point with reference to image blocks;

S5.4 sends static map services request and downloads corresponding to image blocks;

The reference image blocks resampling of downloading is the resolution close with image unit to be matched by S5.5.

Wherein, step S6 comprises further:

S6.1 carries out SIFT coupling to image unit to be matched with reference to image blocks;

S6.2 is by dimensional constraints excluding gross error point;

S6.3 utilizes the anglec of rotation to retrain excluding gross error;

S6.4 utilizes RANSAC to estimate similarity transformation constraint excluding gross error;

S6.5 utilizes affined transformation to retrain excluding gross error point.

Accompanying drawing explanation

Fig. 1 is the on-line automatic matching process process flow diagram according to one embodiment of the present invention;

Fig. 2 is imagery zone and image unit schematic diagram in the inventive method;

Fig. 3 obtains network video map reference point image blocks process flow diagram in the on-line automatic matching process according to one embodiment of the present invention;

Fig. 4 is to image unit to be matched with reference to image blocks coupling process flow diagram in the on-line automatic matching process according to one embodiment of the present invention;

Embodiment

The on-line automatic matching process that the present invention proposes, is described with reference to the accompanying drawings as follows.

As shown in Figure 1, the automatic matching method according to one embodiment of the present invention comprises step:

Remote sensing image to be matched is evenly divided into several regions by the number at the reference mark S1. gathered as required, and each imagery zone is labeled as untreated state, and Fig. 2 is shown in by imagery zone schematic diagram;

If S2. all imagery zones are all processed, then complete and terminate whole Image Matching process; Otherwise, start to process next imagery zone;

S3. current image region is divided into several image units by the size of 256 pixel × 256 pixels, and each image unit is designated as untreated state, Fig. 2 is shown in by image unit schematic diagram;

If all image units S4. in current image region are all processed, then mark this imagery zone and complete process, and proceed to step S2; Otherwise, start to process next image unit;

S5. according to the approximate range with reference to image blocks in the scope of current image unit to be matched and initial imaging model computational grid photomap (network video map available at present comprises Google satellite image map, Bing Aerial image map, MapQuest satellite photomap and Mapbox satellite photomap), then by reference image blocks corresponding to web download and be the resolution close with image unit to be matched by its resampling;

S6. utilize SIFT to mate operator to mate to image unit to be matched with reference to image blocks, and excluding gross error point, if this step obtains the match point of more than 4, then carry out step S7, otherwise proceed to step S4;

S7. from the match point that step S6 obtains, choose a pair that local contrast is maximum, and SIFT is mated the local geometric transformation model obtained as initial value, to match point, Least squares matching is carried out to this, the position coordinate of SIFT feature of refining point, finally this is added in matching result to match point, Least squares matching represents by formula (1) certain any equation of condition

k ₁I _s(x _s，y _s)+k ₂-I _r(x _r，y _r)＝0(1)

Wherein x _s, y _sfor image unit pixel coordinate to be matched, x _r, y _rfor reference image blocks pixel coordinate,

x _r＝a ₀+a ₁x _s+a ₂y _s，y _r＝b ₀+b ₁x _s+b ₂y _s，

A ₀, a ₁, a ₂, b ₀, b ₁, b ₂be 6 geometric transformation parameters, k ₁, k ₂be 2 radiation conversion parameters,

I _s(x _s, y _s) and I _r(x _r, y _r) be respectively image unit to be matched and the gray-scale value with reference to image blocks,

Treat each point in the region of 11 pixel × 11 pixels around match point and set up error equation according to formula (1), then optimization is carried out with Levenberg-Marquardt algorithm, obtain optimum geometric transformation parameter, thus determine the match point coordinate after refining.

After image to be matched is divided into image unit, need according to the approximate range with reference to image blocks in the scope of current image unit to be matched and initial imaging model computational grid photomap.Particularly, as shown in Figure 3, step S5 comprises further:

S5.1 closes on level of zoom most according to the resolution computational grid photomap of image unit to be matched, level of zoom utilizes formula (2) to calculate according to the latitude and longitude coordinates at image unit place to be matched place and resolution, the longitude λ at image unit place place and latitude φ is calculated by the initial imaging model of image to be matched

$\begin{array}{c}GSD=\frac{2{\mathrm{\πR}}_{earth}}{512\×2^{n-1}}\cos \mathrm{\φ}\\ n=\[{\log }_2\frac{2{\mathrm{\πR}}_{earth}\cos \mathrm{\φ}}{512\×GSD}+1\]\end{array}---\left(2\right)$

R in formula _earthrice is earth radius, and its approximate value is 6378137 meters,

GSD is the resolution at image unit place place,

N is level of zoom,

[] represents the computing of getting the most contiguous integer;

S5.2 calculates the corresponding width width with reference to image blocks and height height: given level of zoom n, network video map image coordinate x, the conversion method of y and longitude λ, latitude φ such as formula shown in (3) and formula (4),

$\left\{\begin{array}{c}x=\frac{\mathrm{\λ}+180}{360}\×512\×2^{n-1}\\ y=\left(0.5-\frac{1}{4\mathrm{\π}}\ln \frac{1-\sin \mathrm{\φ}}{1+\sin \mathrm{\φ}}\right)\×512\×2^{n-1}\end{array}\right.---\left(3\right)$

$\left\{\begin{array}{c}\mathrm{\λ}=\frac{360}{512\×2^{n-1}}x-180\\ \mathrm{\φ}=\frac{180}{\mathrm{\π}}arcsin\frac{\exp \[4\mathrm{\π}(0.5-\frac{y}{512\×2^{n-1}})\]-1}{\exp \[4\mathrm{\π}(0.5-\frac{y}{512\×2^{n-1}})\]+1}\end{array}\right.---\left(4\right)$

The latitude and longitude coordinates through type (3) on image unit 4 summits to be matched is utilized to calculate corresponding 4 picture point coordinates respectively, and then ask the minimum enclosed rectangle of these 4 points, the width of this rectangle and the width width be highly with reference to image blocks and height height;

S5.3 calculates the latitude and longitude coordinates λ of the corresponding central point with reference to image blocks _rc, φ _rc, calculate its center dot image coordinate according to the minimum enclosed rectangle obtained in S5.2, then substitute into formula (4) and calculate corresponding latitude and longitude coordinates λ _rc, φ _rc;

S5.4 sends static map services request and downloads corresponding reference image blocks, and static map services request sends in the mode of URL(uniform resource locator) (URL), and the static map services request form of conventional network video map is as follows:

(1) Google satellite image map:

https：//maps.***apis.com/maps/api/staticmap？maptype＝satellite&zoom＝{zoomLevel}&center＝{lat}，{lon}&size＝{width}x{height}&key＝{***Key}

(2) Bing Aerial image map:

http：//dev.virtualearth.net/REST/v1/Imagery/Map/Aerial/{lat}，{lon}/{zoomLevel}？mapSize＝{width}，{height}&key＝{BingMapsKey}

(3) MapQuest satellite photomap:

http：//www.mapquestapi.com/staticmap/v4/getmap？type＝sat&zoom＝{zoomLevel}&center＝{lat}，{lon}&size＝{width}，{height}&key＝{mapquestKey}

(4) Mapbox satellite photomap:

http：//api.tiles.mapbox.com/v4/mapbox.satellite/{lon}，{lat}，{zoomLevel}/{width}x{height}.png？access_token＝{mapboxKey}

Above in each URL, the content in " { } " is filled according to the result of calculation of step S5.1, S5.2 and S5.3, particularly,

ZoomLevel be calculate in step S5.1 close on level of zoom n most,

Width, height are width width and the height height of the reference image blocks calculated in step S5.2,

Lon, lat are the central point longitude λ of the reference image blocks calculated in step S5.3 _rcwith latitude φ _rc,

GoogleKey, BingMapsKey, mapquestKey, mapboxKey are respectively the api key of each Map Service of Network, can apply for obtaining on corresponding website;

The reference image blocks resampling of downloading is the resolution close with image unit to be matched by S5.5.

Download obtains corresponding with reference to after image blocks, needs to utilize SIFT to mate operator and mates to image unit to be matched with reference to image blocks, and excluding gross error point.Particularly, as shown in Figure 4, step S6 comprises further:

S6.1 is respectively from image unit to be matched with reference to extracting SIFT feature point image blocks and calculating corresponding SIFT description vectors, distance metric using the Euclidean distance of SIFT description vectors as two SIFT feature points, two conditions are below judged to each SIFT feature point in image unit to be matched, meet the two one of or meet 2 conditions simultaneously and namely it can be used as matching candidate point:

(1) this is less than 0.75 to the ratio with reference to each SIFT feature point minor increment in image blocks and time small distance;

(2) distance of this point is all not less than with this point apart from the distance of other SIFT feature points in minimum unique point to image unit to be matched with reference in image blocks;

The match point of its correspondence is with reference to SIFT feature point nearest with it in image blocks;

S6.2 checks the SIFT yardstick coordinate ratio T of each candidate point and match point thereof _σif, 0.8≤T _σ≤ 1.25, then retain this to candidate point, otherwise reject this to candidate point;

S6.3 calculates the SIFT principal direction differential seat angle of each candidate point and match point thereof, with 10 degree for setting up grey level histogram in interval, using the angle at histogram peak place as image unit to be matched with reference to the rotation angle θ between image blocks _peak, then check the SIFT principal direction differential seat angle Δ θ of each candidate point and match point thereof, if | Δ θ-θ _peak|≤15 °, then retain this to candidate point, otherwise reject this to candidate point;

S6.4 utilizes RANSAC algorithm to estimate to remain the similarity transformation relation (5) between candidate matches point, rejects the rough error point not meeting this variation relation simultaneously,

$\left\{\begin{array}{c}{x}_r=s\left(x_s\cos \mathrm{\θ}+y_s\sin \mathrm{\θ}\right)+t_x\\ y_r=s\left(-x_s\sin \mathrm{\θ}+y_s\cos \mathrm{\θ}\right)+t_y\end{array}\right.---\left(5\right)$

Wherein x _s, y _sfor image unit pixel coordinate to be matched, x _r, y _rfor reference image blocks pixel coordinate, s and θ is yardstick and the rotation angle parameter of similarity transformation, t _x, t _ybe respectively the translation parameters of similarity transformation in x direction and y direction;

S6.5 calculates affine Transform Model to candidate point remaining in step S6.4,

$\left\{\begin{array}{c}{x}_r=a_0+a_1{x}_s+a_2{y}_s\\ y_r=b_0+b_1{x}_s+b_2{y}_s\end{array}\right.---\left(6\right)$

Wherein x _s, y _sfor image unit pixel coordinate to be matched, x _r, y _rfor reference image blocks pixel coordinate, a ₀, a ₁, a ₂, b ₀, b ₁, b ₂for 6 parameters of affined transformation;

Utilize the affine Transform Model obtained to check each pair of matching candidate point, if maximum residul difference is greater than 1 pixel, then rejects this to after match point, recalculate affine Transform Model, and residual test is carried out to residue candidate point; If maximum residul difference is not more than 1 pixel, then export remaining matching double points, end step S6.

Claims (4)

1., for an on-line automatic matching process for remotely sensing image geometric correction, it is characterized in that: comprise the following steps,

Remote sensing image to be matched is evenly divided into several regions by the number at the reference mark that step S1. gathers as required;

If all imagery zones of step S2. are all processed, then complete whole Image Matching process; Otherwise, start to process next imagery zone;

Current image region is divided into several image units by the size of 256 pixel × 256 pixels by step S3.;

If all image units in step S4. current image region are all processed, then mark this imagery zone and complete process, and proceed to step S2; Otherwise, start to process next image unit;

Step S5: according to the approximate range with reference to image blocks in the scope of current image unit to be matched and initial imaging model computational grid photomap, then by reference image blocks corresponding to web download and be the resolution close with image unit to be matched by its resampling;

Step S6: utilize SIFT to mate operator and mate to image unit to be matched with reference to image blocks, and excluding gross error point, if this step obtains the match point of more than 4, then carry out step S7, otherwise proceed to step S4;

Step S7. selects optimal matching points and carries out Least squares matching to it from the match point that step S6 obtains, and revises the position coordinate of SIFT feature point, this is added in matching result match point.

2. On-line matching method as claimed in claim 1, it is characterized in that, step S5 comprises further:

Step S5.1 closes on level of zoom most according to the resolution computational grid photomap of image unit to be matched;

Step S5.2 calculates the corresponding width with reference to image blocks and height;

Step S5.3 calculates the latitude and longitude coordinates of the corresponding central point with reference to image blocks;

Step S5.4 sends static map services request and downloads corresponding to image blocks;

The reference image blocks resampling of downloading is the resolution close with image unit to be matched by step S5.5.

3. On-line matching method as claimed in claim 1, it is characterized in that, step S6 comprises further:

S6.1 carries out SIFT coupling to image unit to be matched with reference to image blocks;

S6.2 is by dimensional constraints excluding gross error point;

S6.3 utilizes the anglec of rotation to retrain excluding gross error;

S6.4 utilizes RANSAC to estimate similarity transformation constraint excluding gross error;

S6.5 utilizes affined transformation to retrain excluding gross error point.

4. the step S6 of On-line matching method as claimed in claim 3, it is characterized in that, in step S6.1, respectively from image unit to be matched with reference to extracting SIFT feature point image blocks and calculating corresponding SIFT description vectors, distance metric using the Euclidean distance of SIFT description vectors as two SIFT feature points, two conditions are below judged to each SIFT feature point in image unit to be matched, meet the two one of or meet 2 conditions simultaneously and namely it can be used as matching candidate point:

1. this is less than 0.75 to the ratio with reference to each SIFT feature point minor increment in image blocks and time small distance;

2. the distance of this point is all not less than with this point apart from the distance of other SIFT feature points in minimum unique point to image unit to be matched with reference in image blocks.