CN103077527A - Robust multi-source satellite remote sensing image registration method - Google Patents

Abstract

The invention discloses a robust multi-source satellite remote sensing image registration method based on mutual information and block random sample consensus. The method includes the following steps: (A) each layer of pyramid images are generated, and feature points are extracted; (B) on the highest layer of pyramid images, the global registration method is utilized to calculate affine transformation coefficients between the reference image and the slave image, and a rotation angle and a resolution difference coefficient between the images are estimated; (C) the initial positions of homonymy points are accurately predicted, geometric rough correction is carried out on matching window images, and normalized mutual information is utilized to search homonymy points; (D) a quadratic polynomial and a block RANSAC algorithm are utilized to eliminate false matching points; (E) step C and step D are repeated until the original image layer, and accurate image registration is implemented on the basis of the linear rubber sheeting method. The method greatly reduces the workload of manual editing in homonymy point measurement, increases the multi-source satellite remote sensing image data processing automation degree, and can bring remarkable economic and social benefits.

Description

A kind of sane multi-source Satellite Remote Sensing Image method for registering

Technical field

The invention belongs to the Photogrammetry and Remote Sensing technical field, relate to a kind of sane multi-source Remote Sensing Images method for registering.The method can effectively solve the difficult problem that the multi-source Satellite Remote Sensing Image same place is difficult to automatic acquisition based on mutual information Image Matching and piecemeal stochastic sampling consistance technology, realizes the accuracy registration of multi-source Satellite Remote Sensing Image.

Technical background

Along with the develop rapidly of earth observation technology, satellite sensor is more and more diversified.In recent years, High Resolution Remote Sensing Satellites continues to bring out, and Optical remote satellite has Ikonos, Quickbird, and GeoEye and WorldView etc., the SAR remote sensing satellite has TerraSAR-X, Cosmo-Skymed and Radarsat-2 etc.The high-resolution satellite sensor has been widely used in the fields such as topographic mapping, resource exploration, Disaster Assessment.Utilize multiple sensors difference constantly, same type or different types of remote sensing image of the areal that obtains of diverse location carry out the extraction of spatial information, have broad application prospects.

The multi-source Satellite Remote Sensing Image accuracy registration is included it in unified coordinate basis in, is the crucial prerequisite of many application.At present, for the registration of multi-source Remote Sensing Images, the methods that adopt artificial measurement same place more, not only consuming time but also require great effort, can't satisfy again the demand that a large amount of points measure.

For the autoregistration of multi-source Satellite Remote Sensing Image, obtained larger progress.People have proposed the whole registration algorithm, and it takes full advantage of all image informations of image subject to registration, steadily and surely estimate geometric transformation parameter between image by optimization algorithm, and obtained immense success in the medical image registration.The method reliable results, shortcoming are that calculated amount is larger, are only applicable to the registration between small size image.Local registration Algorithm is then mainly chosen the feature of the same name of image to be matched, such as same place, line or face feature, carries out Image registration based on these features of the same name.O.THEPAUT utilizes orbit information that SPOT and ERS image are carried out geometric correction, then extracts respectively edge feature on two width of cloth images, by can realize the registration of two width of cloth images to its coupling.M.D.PAUL proposes a kind of optical image and SAR image matching method based on provincial characteristics, utilizes the description provincial characteristicss such as area, girth, radius and implements the coupling of image.A kind of multi-source multidate high resolution ratio satellite remote-sensing image automatic matching method (number of patent application: 20120296081.5) that Ji Shunping etc. propose, the gradient of at first calculating image blocks obtains gradient image, to the fringe region in the image blocks gradient image and non-fringe region weighting, the non-linear brightness of image blocks gradient image is relevant after the calculating weighting, obtains match point.But, coupling for SAR image and optical image, because both have larger radiation difference, the edge feature that extracts on the optical image is abundant and continuous, but the SAR image is because the impact of speckle noise is difficult to obtain complete edge, this has brought very large difficulty for the coupling of feature of the same name, and the method for M.D.PAUL then requires image to have obvious planar feature, lacks broad applicability.Zhang Hong etc. have proposed high resolution SAR image registration processing method and have applied for patent that (number of patent application: 200510132200.3), but it is mainly for homology SAR Image registration, does not consider that multi-source Remote Sensing Images radiation significant difference is difficult to the problems such as coupling.

Stochastic sampling consistance (Random Sample Consensus, RAN SAC), it can estimate from the data centralization that comprises a large amount of erroneous point high-precision parameter, this patent makes improvements, use it for detection and the rejecting of error matching points in the large format multi-source Remote Sensing Images coupling, and merge mutual information coupling and pyramid matching strategy, finally realize the accuracy registration of multi-source Satellite Remote Sensing Image.

Summary of the invention

In present multi-source Satellite Remote Sensing Image autoregistration process, exist error matching points more, reject the Mismatching point shortcoming such as waste time and energy manually, the objective of the invention is to be to provide a kind of based on mutual information and the conforming multi-source Remote Sensing Images method for registering of piecemeal stochastic sampling, realize automatic, reliable, the exact matching of same place, greatly improve the automaticity of satellite image registration, greatly reduce labor workload.

In order to achieve the above object, the present invention adopts following technical measures:

A kind of sane multi-source Satellite Remote Sensing Image method for registering the steps include:

A, utilize 3 * 3 pixel methods of average, generate 3 layers of pyramid image, and utilize Forstner feature extraction operator, 30 * 30 graticule mesh that evenly distribute on main image subject to registration are extracted a unique point from each graticule mesh.

B, at top pyramid image, take the normalized mutual information of major and minor image subject to registration as the target equation, seek the affined transformation coefficient vector of global registration by Powell optimization (books " Numerical recipes in C " that can publish with reference to Cambridge University Press 1992), on this basis, calculate the anglec of rotation and differences in resolution coefficient approximate between image.

Before C, the Feature Points Matching, at first need Accurate Prediction same place initial position.When top pyramid image coupling, directly utilize the affined transformation coefficient that calculates among the step B to carry out a position prediction; When other pyramid image couplings, then utilize upper strata pyramid image matching result, the anglec of rotation between the consideration image and the impact of differences in resolution according to the continuous characteristic of neighbor point parallax, are carried out the prediction of same place initial position.Then, the match window image is carried out the thick correction of geometry, central point corresponding to normalized mutual information maximized window is same place in the search window.

D, the view picture satellite remote-sensing image is divided into some, such as 3 * 3,5 * 5, the geometrical-restriction relation of matching result is expressed with quadratic polynomial in each piece, utilizes successively RANSAC algorithm (can with reference to the article " Random Sample Consensus:A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography " of M.A.Fischler and R.C.Bolles) to delete error matching points in every of this layer pyramid image matching result.

E, repeating step C and D are until finish the raw video layer, utilize at last the major and minor image same place that obtains, adopt linear rubber pulling method, image is divided into some little triangles, adopt the affined transformation method to resample to each triangle, finally realize the accuracy registration of satellite remote-sensing image.

Multi-group data is tested, first group is two scape Envisat SAR images, second group is COSMO-SkyMed SAR image and IRS-P5 image, all can automatically steadily and surely match more equally distributed unique point, realize the accurate Pixel-level registration of satellite remote-sensing image, line feature of the same name etc. all can well be coincide on the image.

The present invention combines the methods such as match window image geometric rough correct, normalized mutual information coupling, piecemeal RANSAC algorithm, pyramid image coupling, can realize automatically reliably obtaining of multi-source Satellite Remote Sensing Image same place, and then finish its accuracy registration.The method has greatly reduced human-edited's amount that same place measures, and improves the datamation degree of multi-source Satellite Remote Sensing Image, has significant economic benefit and social benefit.

Compared with prior art, remarkable advantage of the present invention and effect are mainly manifested in:

(1) during successively pyramid image mates, can automatically find and reject error matching points in each layer matching result, have better coupling reliability;

(2) whether satellite remote-sensing image is evenly distributed obvious wire, planar feature without specific (special) requirements, have stronger adaptability.

The present invention provides a kind of reliable same place matching process for the automatic accuracy registration of multi-source Satellite Remote Sensing Image, introduce piecemeal RANSAC algorithm in pyramid image coupling successively, can effectively delete the error matching points in the matching result, realize automatic, reliable, the exact matching of same place, improve the automaticity of multi-source Satellite Remote Sensing Image registration.

Description of drawings

Fig. 1 is a kind of process flow diagram of the multi-source Satellite Remote Sensing Image method for registering based on mutual information and piecemeal RANSAC algorithm.

Embodiment

Below in conjunction with accompanying drawing the present invention is done and to describe in further detail.

Embodiment 1:

A kind of multi-source Satellite Remote Sensing Image method for registering based on mutual information and piecemeal RANSAC algorithm, its process flow diagram as shown in Figure 1, each step is elaborated as follows:

The first step, the data pre-service 1 of multi-source Satellite Remote Sensing Image coupling:

The data pre-service of multi-source Satellite Remote Sensing Image coupling mainly comprises: the generation of pyramid image, the extraction of unique point.

Concrete steps are as follows:

The generation of 1 pyramid image: adopt practical, simple 3 * 3 pixel methods of average, will survey 3 grades of pyramid images of all video generations of district.At first per 3 * 3 pixels of raw video (0 grade of pyramid image) calculated its average gray value, and be assigned to the corresponding pixel of the 1st grade of pyramid image, generate first order pyramid image.The rest may be inferred, until generate the 2nd grade of pyramid image.

The extraction of 2 unique points: image is divided into uniform graticule mesh, graticule mesh is counted can be given according to the image size, generally image is divided into 30 * 30 and gets final product with interior graticule mesh, adopt photogrammetric in the feature extraction operators such as Forstner commonly used in each graticule mesh, extract a best features point.If the characteristic information in certain graticule mesh is not obvious, then with the graticule mesh central point as unique point.Like this, can guarantee that not only unique point is evenly distributed on whole image, and can set up easily the proximity relations between the unique point.

Second step, the prediction 2 of unique point initial point position:

In order to reduce the search time of same place, improve the success ratio of Image Matching, need to carry out Accurate Prediction to unique point initial point position.

Concrete steps are as follows:

The same place prediction of 1 top pyramid image coupling: on top pyramid image, at first, by the global registration algorithm, utilize optimization algorithm to find the solution geometric transformation parameter between image subject to registration.On top pyramid image, it is milder that parallax changes, and can utilize affined transformation to come geometric relationship between the approximate expression image.

Specific as follows:

x′＝p ₁+p ₂x+p ₃y

(1)

y＝p ₄+p ₅x+p ₆y

(x ', y '), (x, y) represent respectively the picpointed coordinate of same place on main image (reference image) and subpictures (slave image), p in the formula ₁, p ₂, p ₃, p ₄, p ₅, p ₆Be affine transformation parameter to be found the solution.

The global registration algorithm is in fact the process of an iteration optimizing, when target equation F value maximum, the at this moment parameter of gained

Be required, namely

$\hat{p}=\arg \underset{p}{\max }F(R,S;p)---\left(2\right)$

In the formula, F is the target equation, and concrete meaning is the normalization information (Normalized Mutual Information, NMI) of main image and subpictures in the inventive example, and R is main image, and S is subpictures, and p is affined transformation coefficient vector (p to be found the solution ₁, p ₂, p ₃, p ₄, p ₅, p ₆) ^T

Expression makes the searching process of the parameter p of target equation F maximum,

Expression makes the affined transformation coefficient vector of major and minor image normalized mutual information maximum.

The optimized algorithm of parameter search has a lot, selected the Powell optimization among the present invention, the method does not need objective function is carried out derived function, has advantages of direct method, and have quadratic convergence, be acknowledged as the highly effective method of present solution Unconstrained Optimization Problem.

The computing formula of normalized mutual information NMI is as follows:

$\mathrm{NMI}(A,B)=\frac{H\left(A\right)+H\left(B\right)}{H(A,B)}---\left(3\right)$

In the formula, A, B represent respectively image blocks to be matched, and H (A) and H (B) are respectively the Averages of image block A and B, and H (A, B) is the joint entropy of two image blocks.

The Average Mutual of A and B and joint entropy can followingly be calculated:

$H\left(A\right)=\underset{a}{\mathrm{\Σ}}-P_A\left(a\right)\log P_A\left(a\right)$

$H\left(B\right)=\underset{b}{\mathrm{\Σ}}-P_B\left(b\right)\log P_B\left(b\right)---\left(4\right)$

$H(A,B)=\underset{a,b}{\mathrm{\Σ}}-P_{A,B}(a,b)\log P_{A,B}(a,b)$

In the formula, P _A(a) and P _B(b) be marginal probability density function, P _{A, B}(a, b) is joint probability density function.They can be obtained by following formula:

$P_A\left(a\right)=\underset{a}{\mathrm{\Σ}}{P}_{A,B}(a,b)$

$P_B\left(b\right)=\underset{b}{\mathrm{\Σ}}{P}_{A,B}(a,b)---\left(5\right)$

$P_{A,B}(a,b)=\frac{h(a,b)}{\underset{a,b}{\mathrm{\Σ}}h(a,b)}$

In the formula, h represents the joint histogram of image blocks A and B.It is a two-dimensional matrix:

Here, h (a, b) is the element of joint histogram h, has gray-scale value a among its presentation video piece A, has the point of gray-scale value b to number in image block B.M and N be the number of greyscale levels of presentation video piece A, image block B respectively.In the inventive example, in order to reduce calculated amount, the original gray scale of image is stretched, M and N all get 32.

Calculate the affined transformation coefficient, utilize affine Transform Model can obtain the initial position of same place on subpictures of main image feature point in the top pyramid.

The same place prediction of 2 other layers pyramid image coupling: at other pyramid image layers, then utilize the upper strata matching result, according to the continuous characteristic of neighbor point parallax, the same place initial position is predicted.

Concrete steps are as follows:

1) anglec of rotation and the differences in resolution coefficient between the calculating image: utilize 6 affined transformation coefficients that the global registration algorithm calculates in the top pyramid image coupling, the anglec of rotation between image and differences in resolution coefficient are carried out approximate evaluation.

For the anglec of rotation between image, at first calculate two valuation α of the anglec of rotation ₁, α ₂:

$\left\{\begin{array}{c}{\mathrm{\α}}_1=\arctan (a_2/a_1)\\ {\mathrm{\α}}_2=\arctan (-b_1/b_2)\end{array}\right.---\left(7\right)$

In the formula, arctan represents the arctan function computing in the trigonometric function.With α ₁, α ₂Mean value can be used as the valuation of the two image anglecs of rotation

In like manner, at first calculate two valuation λ of the differences in resolution coefficient between image ₁, λ ₂:

$\left\{\begin{array}{c}{\mathrm{\λ}}_1=\sqrt{a_1^2+a_2^2}\\ {\mathrm{\λ}}_2=\sqrt{a_3^2+a_4^2}\end{array}\right.---\left(8\right)$

With λ ₁, λ ₂Mean value can be used as the valuation of two width of cloth image resolution coefficients of variation

2) for the upper strata point that the match is successful, directly projection obtains the initial position at this layer; But for the point that it fails to match, because it utilizes neighbor point parallax future position position, the impact that the anglec of rotation between needs consideration image and differences in resolution are brought.

Suppose p ₁(x ₁, y ₁) be that it fails to match on main image upper strata point, p ₂(x ₂, y ₂) be it is contiguous the match is successful point, p ' ₂(x ' ₂, y ' ₂) be p ₂Identical point coordinates on subpictures.

At first calculate p ₁And p ₂2 at x, the displacement on the y direction

$\left\{\begin{array}{c}\mathrm{dx}=x_2-x_1\\ \mathrm{dy}=y_2-y_1\end{array}\right.---\left(9\right)$

Then calculate the displacement size behind the rotation alpha angle

$\left\{\begin{array}{c}{\mathrm{dx}}^{\′}=\mathrm{dx}\cos \mathrm{\α}-\mathrm{dy}\sin \mathrm{\α}\\ {\mathrm{dy}}^{\′}=\mathrm{dx}\sin \mathrm{\α}+\mathrm{dy}\cos \mathrm{\α}\end{array}\right.---\left(10\right)$

Displacement after the calculating compensation differences in resolution is big or small again

$\left\{\begin{array}{c}{\mathrm{dx}}^{\′\′}=\mathrm{\λ}{\mathrm{dx}}^{\′}\\ {\mathrm{dy}}^{\′\′}=\mathrm{\λ}{\mathrm{dy}}^{\′}\end{array}\right.---\left(11\right)$

P then ₁Initial position (the x of same place on subpictures of point ₁', y ₁') be:

$\left\{\begin{array}{c}{x}_1^{\′}=x_2^{\′}+{\mathrm{dx}}^{\′\′}\\ y_1^{\′}=y_2^{\′}+{\mathrm{dy}}^{\′\′}\end{array}\right.---\left(12\right)$

The 3rd step, how much thick normalized mutual information couplings 3 of correcting of band

In order to reduce the adverse effect of the image anglec of rotation and geometry deformation, utilize normalized mutual information to estimate before the search same place, the match window image is carried out the thick correction of geometry.

Concrete steps are as follows:

The geometry of 3 match window images is slightly corrected: in the inventive example, the search window image of directly treating matching characteristic point carries out the thick correction of geometry, when the pixel in the traversal search window is mated like this, need not to resample to obtain the match window image at every turn.

Take the initial position A of unique point to be matched on subpictures as initial point, set up local coordinate system A-x ' y ' and A-xy, wherein x axle, y axle are that the respective coordinates axle is parallel to each other with the subpictures pixel coordinate respectively, x ' axle, y ' axle respectively with rotation compensation after search window respective coordinates axle be parallel to each other, then the corner dimension between two coordinate systems is approximately the anglec of rotation α between image.Then, utilize formula (13) can calculate the coordinate of pixel in local coordinate system A-xy in the search window:

$\left\{\begin{array}{c}x=(x^{\′}\cos \mathrm{\α}-y^{\′}\sin \mathrm{\α})\×\mathrm{\λ}\\ y=(x^{\′}\sin \mathrm{\α}+y^{\′}\cos \mathrm{\α})\×\mathrm{\λ}\end{array}\right.---\left(13\right)$

Then, add pixel coordinate corresponding to initial point position A, can obtain the pixel coordinate of this pixel on subpictures in the search window.By bilinear interpolation method, can calculate corresponding gray-scale value.Each pixel in the traversal processing search window, the final generation through the search window image after the anglec of rotation and the differences in resolution compensation.

4 normalized mutual informations coupling: the geometry that carries out the match window image by upper one step process is slightly corrected, and then utilizes formula (3) to calculate normalized mutual information.In the inventive example, the match window size is 61 * 61 pixels.At last, choose the pixel of normalized mutual information value maximum in the search window as the same place of this unique point.

In the 4th step, reject error matching points 4 based on piecemeal RANSAC:

In general, remote sensing satellite operates on the track of hundreds of kilometer height, and the image point displacement less that topographic relief causes exists stronger systematicness, common available quadratic polynomial model description between the picpointed coordinate of image to be matched.But because the breadth of satellite remote-sensing image is often larger, and for high resolving power mountain area image, its parallax changes relatively violent, and whole match point of rejecting mistake with a quadratic polynomial model tends to correct matching result mistake is deleted.

The piecemeal RANSAC that the present invention proposes rejects the error matching points method, at first the view picture remote sensing image is divided into some, such as 3 * 3 piecemeals or 4 * 4 piecemeals or 5 * 5 piecemeals, then the match point in every is used respectively its geometrical-restriction relation of quadratic polynomial model description, and auxiliary come automatic detecting and elimination error matching points in the RANSAC algorithm.

Quadratic polynomial model concrete form is:

$\left\{\begin{array}{c}{x}_s=a_0+a_1{x}_m+a_2{y}_m+a_3{y}_m+a_4{x}_m^2+a_5{y}_m^2\\ y_s=b_0+b_1{x}_m+b_2{y}_m+b_3{x}_m{y}_m+b_4{x}_m^2+b_5{y}_m^2\end{array}\right.---\left(14\right)$

In the formula, (x _m, y _m), (x _s, y _s) be respectively the pixel coordinate of corresponding image points on main image and subpictures.The solution procedure of RANSAC algorithm is as follows:

1) establish sample number k for infinitely great, Sample Counter t is 0;

2) from current piecemeal, randomly draw 6 pairs of points in the matching characteristic point and set up system of equations, resolve 12 coefficients of quadratic polynomial;

3) calculate in the current piecemeal all the other unique points through behind the polynomial transformation and the distance between its match point.If distance is thought that then this correctly mates same place, otherwise is erroneous matching less than given threshold value.The last shared ratio of mistake of statistics matching double points

4) calculate sample numerical value

η gets 0.99 in the inventive example;

5) Sample Counter t adds 1;

6) when k＜t, stopping circulation, otherwise rotate back into the 2nd) step recomputates;

7) after iteration stops, choose the set with at most correct matching double points and recomputate the quadratic polynomial coefficient, and add up each to the error of fitting of same place, disallowable greater than being considered as error matching points of given threshold value.

In the 5th step, carry out image accuracy registration 5 based on linear rubber pulling method

At first with the same place on all major and minor images to forming the triangulation network, to each delta-shaped region, calculate the affined transformation coefficient by least-squares algorithm, then each delta-shaped region utilizes, utilize corresponding affined transformation coefficient to carry out coordinate conversion, calculate gray-scale value by bilinear interpolation, finally finish the accuracy registration between image.

Claims (1)

1. a sane multi-source Satellite Remote Sensing Image method for registering the steps include:

The data pre-service (1) of A, multi-source Satellite Remote Sensing Image coupling:

The data pre-service of multi-source Satellite Remote Sensing Image coupling comprises: the generation of pyramid image, the extraction of unique point: step is:

The generation of a, pyramid image: adopt 3 * 3 pixel methods of average, to survey 3 grades of pyramid images of all video generations of district, at first per 3 * 3 pixels of raw video calculated its average gray value, and be assigned to the corresponding pixel of the 1st grade of pyramid image, generate first order pyramid image, until generate the 2nd grade of pyramid image;

The extraction of b, unique point: image is divided into uniform graticule mesh, graticule mesh is counted given according to the image size, image is divided into 30 * 30 with interior graticule mesh, adopt photogrammetric in Forstner feature extraction operator commonly used in each graticule mesh, extract a unique point, characteristic information in the graticule mesh is not obvious, with the graticule mesh central point as unique point;

The prediction (2) of B, unique point initial point position:

In order to reduce the search time of same place, improve the success ratio of Image Matching, Accurate Prediction is carried out in unique point initial point position:

Step is as follows:

The same place prediction of a, top pyramid image coupling: on top pyramid image, at first, by the global registration algorithm, geometric transformation parameter between utilizationization Algorithm for Solving image subject to registration, on top pyramid image, it is milder that parallax changes, and utilizes affined transformation to come geometric relationship between the approximate expression image, specific as follows:

x′＝p ₁+p ₂x+p ₃y (1)

y′＝p ₄+p ₅x+p ₆y

(x ', y '), (x, y) represent respectively the picpointed coordinate of same place on main image and subpictures, p in the formula ₁, p ₂, p ₃, p ₄, p ₅, p ₆Be the affine transformation parameter of finding the solution;

The global registration algorithm is the process of an iteration optimizing, and target equation F value is maximum, the parameter of gained Be required, namely

$\hat{p}=\arg \underset{p}{\max }F(R,S;p)---\left(2\right)$

In the formula, F is the target equation, and R is main image, and S is subpictures, and p is affined transformation coefficient vector (p to be found the solution ₁, p ₂, p ₃, p ₄, p ₅, p ₆) ^T,

Expression makes the searching process of the parameter p of target equation F maximum,

Expression makes the affined transformation coefficient vector of major and minor image normalized mutual information maximum;

The computing formula of normalized mutual information NMI is as follows:

$\mathrm{NMI}(A,B)=\frac{H\left(A\right)+H\left(B\right)}{H(A,B)}---\left(3\right)$

In the formula, A, B represent respectively image blocks to be matched, and H (A) and H (B) are respectively the Averages of image block A and B, and H (A, B) is the joint entropy of two image blocks;

The following calculating of the Average Mutual of A and B and joint entropy:

$H\left(A\right)=\underset{a}{\mathrm{\Σ}}-P_A\left(a\right)\log P_A\left(a\right)$

$H\left(B\right)=\underset{b}{\mathrm{\Σ}}-P_B\left(b\right)\log P_B\left(b\right)---\left(4\right)$

$H(A,B)=\underset{a,b}{\mathrm{\Σ}}-P_{A,B}(a,b)\log P_{A,B}(a,b)$

In the formula, P _A(a) and P _B(b) be marginal probability density function, P _{A, B}(a, b) is joint probability density function, obtained by following formula:

$P_A\left(a\right)=\underset{a}{\mathrm{\Σ}}{P}_{A,B}(a,b)$

$P_B\left(b\right)=\underset{b}{\mathrm{\Σ}}{P}_{A,B}(a,b)---\left(5\right)$

$P_{A,B}(a,b)=\frac{h(a,b)}{\underset{a,b}{\mathrm{\Σ}}h(a,b)}$

In the formula, h represents the joint histogram of image blocks A and B, is a two-dimensional matrix:

Here, h (a, b) is the element of joint histogram h, has gray-scale value a among the presentation video piece A, the point of gray-scale value b is arranged to number in image block B, and M and N be the number of greyscale levels of presentation video piece A, image block B respectively;

The same place prediction of other layers pyramid image coupling: at other pyramid image layers, utilize the upper strata matching result, according to the continuous characteristic of neighbor point parallax, the same place initial position is predicted;

Step is as follows:

1) anglec of rotation and the differences in resolution coefficient between the calculating image: utilize 6 affined transformation coefficients that the global registration algorithm calculates in the top pyramid image coupling, the anglec of rotation between image and differences in resolution coefficient are carried out approximate evaluation;

For the anglec of rotation between image, at first calculate two valuation α of the anglec of rotation ₁, α ₂:

$\left\{\begin{array}{c}{\mathrm{\α}}_1=\arctan (a_2/a_1)\\ {\mathrm{\α}}_2=\arctan (-b_1/b_2)\end{array}\right.---\left(7\right)$

In the formula, arctan represents the arctan function computing in the trigonometric function, with α ₁, α ₂Mean value as the valuation of the two image anglecs of rotation

In like manner, at first calculate two valuation λ of the differences in resolution coefficient between image ₁, λ ₂:

$\left\{\begin{array}{c}{\mathrm{\λ}}_1=\sqrt{a_1^2+a_2^2}\\ {\mathrm{\λ}}_2=\sqrt{a_3^2+a_4^2}\end{array}\right.---\left(8\right)$

With λ ₁, λ ₂Mean value as the valuation of two width of cloth image resolution coefficients of variation

2) for the upper strata point that the match is successful, directly projection obtains the initial position at this layer; For the point that it fails to match, utilize neighbor point parallax future position position, the impact that the anglec of rotation between the consideration image and differences in resolution are brought;

If p ₁(x ₁, y ₁) be that it fails to match on main image upper strata point, p ₂(x ₂, y ₂) be it is contiguous the match is successful point, p ' ₂(x ' ₂, y ' ₂) be p ₂Identical point coordinates on subpictures;

At first calculate p ₁And p ₂2 at x, the displacement on the y direction

$\left\{\begin{array}{c}\mathrm{dx}=x_2-x_1\\ \mathrm{dy}=y_2-y_1\end{array}\right.---\left(9\right)$

Then calculate the displacement size behind the rotation alpha angle

$\left\{\begin{array}{c}{\mathrm{dx}}^{\′}=\mathrm{dx}\cos \mathrm{\α}-\mathrm{dy}\sin \mathrm{\α}\\ {\mathrm{dy}}^{\′}=\mathrm{dx}\sin \mathrm{\α}+\mathrm{dy}\cos \mathrm{\α}\end{array}\right.---\left(10\right)$

Displacement after the calculating compensation differences in resolution is big or small again

$\left\{\begin{array}{c}{\mathrm{dx}}^{\′\′}=\mathrm{\λ}{\mathrm{dx}}^{\′}\\ {\mathrm{dy}}^{\′\′}=\mathrm{\λ}{\mathrm{dy}}^{\′}\end{array}\right.---\left(11\right)$

p ₁The initial position of same place on subpictures of point (x ' ₁, y ' ₁) be:

$\left\{\begin{array}{c}{x}_1^{\′}=x_2^{\′}+d{x}^{\′\′}\\ y_1^{\′}=y_2^{\′}+d{y}^{\′\′}\end{array}\right.---\left(12\right)$

C, how much thick normalized mutual information couplings (3) of correcting of band:

Reduce the adverse effect of the image anglec of rotation and geometry deformation, utilize normalized mutual information to estimate before the search same place, the match window image is carried out the thick correction of geometry:

Step is as follows:

The geometry of a, match window image is slightly corrected: the search window image of directly treating matching characteristic point carries out the thick correction of geometry, when the pixel in the traversal search window is mated like this, need not to resample to obtain the match window image at every turn;

The initial position A of unique point to be matched on subpictures is initial point, set up local coordinate system A-x ' y ' and A-xy, wherein x axle, y axle are that the respective coordinates axle is parallel to each other with the subpictures pixel coordinate respectively, x ' axle, y ' axle respectively with rotation compensation after search window respective coordinates axle be parallel to each other, corner dimension between two coordinate systems is approximately the anglec of rotation α between image, utilizes formula (13) to calculate the coordinate of pixel in local coordinate system A-xy in the search window:

$\left\{\begin{array}{c}x=(x^{\′}\cos \mathrm{\α}-y^{\′}\sin \mathrm{\α})\×\mathrm{\λ}\\ y=(x^{\′}\sin \mathrm{\α}+y^{\′}\cos \mathrm{\α})\×\mathrm{\λ}\end{array}\right.---\left(13\right)$

Add pixel coordinate corresponding to initial point position A, obtain the pixel coordinate of this pixel on subpictures in the search window, pass through bilinear interpolation method, calculate corresponding gray-scale value, each pixel in the traversal processing search window, the final generation through the search window image after the anglec of rotation and the differences in resolution compensation;

B, normalized mutual information coupling: the geometry that carries out the match window image by upper one step process is slightly corrected, utilize formula (3) to calculate normalized mutual information, the match window size is 61 * 61 pixels, at last, choose the pixel of normalized mutual information value maximum in the search window as the same place of this unique point;

D, reject error matching points (4) based on piecemeal RANSAC:

The remote sensing satellite operation in orbit, the image point displacement that topographic relief causes is relatively little, there is systematicness between the picpointed coordinate of image to be matched, usually use the quadratic polynomial model description, the piecemeal RANSAC that proposes rejects the error matching points method, the view picture remote sensing image is divided into 3 * 3 piecemeals or 4 * 4 piecemeals or 5 * 5 piecemeals, then the match point in every is used respectively its geometrical-restriction relation of quadratic polynomial model description, and auxiliary come automatic detecting and elimination error matching points in the RANSAC algorithm;

Quadratic polynomial model concrete form is:

$\left\{\begin{array}{c}{x}_s=a_0+a_1{x}_m+a_2{y}_m+a_3{x}_m{y}_m+a_4{x}_m^2+a_5{y}_m^2\\ y_s=b_0+b_1{x}_m+b_2{y}_m+b_3{x}_m{y}_m+b_4{x}_m^2+b_5{y}_m^2\end{array}\right.---\left(14\right)$

In the formula, (x _m, y _m), (x _s, y _s) being respectively the pixel coordinate of corresponding image points on main image and subpictures, the solution procedure of RANSAC algorithm is as follows:

1) establish sample number k for infinitely great, Sample Counter t is 0;

2) from current piecemeal, randomly draw 6 pairs of points in the matching characteristic point and set up system of equations, resolve 12 coefficients of quadratic polynomial;

3) calculate in the current piecemeal all the other unique points through behind the polynomial transformation and the distance between its match point, distance is less than given threshold value, and this correctly mates same place, the last shared ratio of mistake of statistics matching double points

4) calculate sample numerical value $k=\frac{\log (1-\mathrm{\η})}{\log (1-{(1-{\mathrm{\ϵ}}^6)}^t)},$ η gets 0.99;

5) Sample Counter t adds 1;

6) during k＜t, stop circulation;

7) after iteration stops, choose the set of at most correct matching double points and recomputate the quadratic polynomial coefficient, and add up each to the error of fitting of same place;

E, carry out image accuracy registration (5) based on linear rubber pulling method:

At first with the same place on all major and minor images to forming the triangulation network, to each delta-shaped region, calculate the affined transformation coefficient by least-squares algorithm, each delta-shaped region utilizes, utilize corresponding affined transformation coefficient to carry out coordinate conversion, calculate gray-scale value by bilinear interpolation, finally finish the accuracy registration between image.

Images