CN101251926B - Remote sensing image registration method based on local configuration covariance matrix - Google Patents

Abstract

The invention provides a method for registering a remote sensing image based on a local contour covariance matrix, which combines a corner feature, a local split image and a local edge contour together as a local feature to implement an extraction of control points, has the local contour covariance matrix introduced into a registering processing, uses the local contour covariance matrix to implement an automatic extraction of a rotating and scaling factor, and combines a technology of self adaptive selection of a local window to reduce an effect of the twiddle and scaling factor on the registering processing so as to improve the robustness of the registering processing, thereby implementing an accurate registering between remote sensing images in a mode of large rotation and scaling without human participation.

Description

A kind of remote sensing image registration method based on local configuration covariance matrix

Technical field

The invention belongs to image processing field, relate to a kind of method for registering images, be specifically related to a kind of remote sensing image registration disposal route based on local configuration covariance matrix.

Background technology

Image registration is a underlying issue of Flame Image Process, is the process that two width of cloth or multiple image to the Same Scene of taking from different time, different sensors or different visual angles mate, superpose.Say that accurately the target of image registration is exactly the best mapping relations that find between the input picture.

At present, image registration disposal route mainly can be divided three classes:

1, based on the registration process method of overall gray-scale statistical information;

The gray-scale statistical characteristic of directly utilizing image itself to have based on the registration process method of overall gray-scale statistical information is come the similarity between the dimensioned plan picture, and then realizes the registration process between image.This disposal route great advantage is that principle is simple, is easy to realize, for obtaining the good treatment effect under the less situation of similar image and affine error.Yet for foreign peoples's remote sensing images, because there is bigger difference in the gradation of image characteristic, its treatment effect is unsatisfactory.Simultaneously, owing to adopt iterative search in the processing procedure, when affine error was big, its processing speed sharply descended, even can't restrain.

2, based on the registration process method of phase place correlation technique;

Theoretical foundation based on the registration process method of phase place correlation technique is the Fourier phase correlation technique, utilizes frequency domain information to carry out relevant treatment, the search optimum matching, and then realize automatic image registration not needing to seek under the situation at reference mark.The advantage of this algorithm is that to insensitive for noise owing to can adopt ready-made fft algorithm, therefore, this registration process algorithm has certain advantage on processing speed.In actual use, this algorithm needs to have bigger overlapping region between image subject to registration, and the gray feature of image is similar, and this has limited the use of this Processing Algorithm in the remote sensing image registration process.

3, based on the registration process method of characteristics of image;

Based on the registration process method of characteristics of image is another big class in the method for registering images, as shown in Figure 1, this method is read in each image from each image channel after, at first extract the characteristic information in the image, utilize the feature of extracting to carry out the local similar coupling, and set up mapping relations between the image by the matching relationship of feature, it is right to handle the image that obtains coupling by image resampling at last.

Comparatively speaking, because the too much gamma characteristic that depends on image of preceding two class disposal routes, their application in remote sensing image registration is handled have been restricted, therefore, registration process method based on feature becomes the common method that remote sensing image registration is handled, and its core is to set up the mapping relations between characteristics of image.

Aspect characteristic matching relation sets up, traditional method be by statistical parameters such as statistics partial auto-correlation or interactive information weigh in the image 2 whether consistent, but this method is subjected to the influence of twiddle factor and zoom factor very big, along with the increase of zoom factor and twiddle factor, its correlativity reduces rapidly.Utilize affine invariant to carry out the similarity test, then depend on the edge contour that extracts complete closure, and this is difficult to usually for remote sensing images especially radar image realize.The point feature of utilizing aggregation between the some feature to realize that the similarity coupling then requires each input picture and extracted has very high consistance, and for remote sensing images, this is difficult to realize equally.Therefore, how effectively to set up the mapping relations between input picture, become the difficult point that remote sensing image registration is handled.

Summary of the invention

The present invention proposes a kind of remote sensing image registration method based on local configuration covariance matrix, with the angle point feature, local segmentation image and local edge contour combine, realize the reference mark extraction as local feature jointly, and local configuration covariance matrix introduced registration process, utilize local configuration covariance matrix to realize rotating the automatic extraction of zoom factor, reduce to rotate the influence of zoom factor in conjunction with local window self-adaptation selecting technology to registration process, improve the robustness of registration process, realizing the accurate registration between remote sensing images under the big rotation zoom mode under the no artificial participation condition.

At first extract the angle point feature in the input picture, be characterized as reference with angle point, extract near the local configuration feature of angle point, and the covariance matrix that utilizes the local configuration feature extracts the rotation zoom factor between input picture, reduce to rotate the influence of zoom factor, improve the extraction precision of same place characteristic matching; Then, utilize detected same place to carry out the affine model parameter estimation; At last, carry out image resampling according to the transformation model parameter that estimates, the remote sensing images behind the acquisition registration are right.On the whole, this treatment scheme automatic extraction, optimum window size of mainly comprising the selection of feature and extraction, rotation zoom factor choose, angle point characteristic matching, affine model parameter estimation and five parts of image resampling.

A kind of remote sensing image registration method based on local configuration covariance matrix comprises the steps:

Step 1: from each image channel, read in each image simultaneously, adopt the Harris Corner Detection Algorithm to realize the image Corner Detection, judge whether to exist angle point by the angular response functional value of adding up each point in the mobile search window;

If the angular response function C RF value of any is the maximal value in the local neighborhood, then this point is set to alternative angle point; Otherwise show that this point is not an angle point.After in having judged image, having a few, utilize threshold processing to extract final angle point.

Step 2: adopt inter-class variance automatic threshold method that reference picture and image subject to registration are carried out local threshold and cut apart; Utilization is extracted to such an extent that optimal threshold is cut apart topography, obtains local bianry image, and obtains the local edge profile according to the local bianry image that extracts.

Step 3: the automatic extraction of rotation zoom factor; Employing is based on the rotation zoom factor extraction method of local configuration covariance matrix.The bianry image that utilizes Threshold Segmentation to obtain calculates the covariance matrix of reference picture and image subject to registration respectively; Distribute according to the partial structurtes of reference picture and image subject to registration and can extract twiddle factor and zoom factor between image.

Step 4: choose the analysis window size by the ratio of selected window size in comparative analysis reference picture and the image subject to registration and the zoom factor that utilizes input picture to extract, extract optimum window.

Step 5: utilize the rotation zoom factor and the optimum window size that extract that the edge contour that extracts in the step 1 is resampled, reduce to rotate the influence of zoom factor to characteristic matching; And the edge feature after resampling carried out the Freeman coding, and add up similarity between two line segments by the chain code correlation technique; Utilize the mutual support characteristic between relative position relation between point set and point set to reject error matching points.

Step 6: after finishing Feature Points Matching, utilize the parameter in the least square method estimation affine model; After determining transformation parameter, according to affine parameter to input picture to processings that resample, it is right to obtain the image that mates mutually.

In the described step 3, at the automatic extraction of local segmentation image stretching rotation zoom factor.It is bright, dark two class images that passing threshold is cut apart image segmentation, adds up two-part covariance matrix and rotation zoom factor then respectively, determines the corresponding relation between image, and extracts corresponding rotation zoom factor.

In the described step 3, when determining the corresponding relation between image, these two kinds of corresponding relations of dark, bright image in bright, dark image in the corresponding image subject to registration of bright in the comparative analysis reference picture, dark image and the reference picture in the corresponding image subject to registration of bright, dark image, choose the ratio of anglec of rotation difference and zoom factor, it is littler to choose anglec of rotation difference, and the more approaching corresponding relation of the ratio of diaxon rotation zoom factor is as the corresponding relation between two images.

In the described step 4, the zoom factor that the statistics local configuration covariance matrix extracts, calculate the ratio that twiddle factor and zoom factor calculate selected window size in comparative analysis reference picture and the image subject to registration then, and choose under the different windows size condition both difference in the comparative analysis image subject to registration, if it is minimum that both differences reach, then obtain the size of the size of this window subject to registration as optimum window; Otherwise recomputate the zoom factor that local configuration covariance matrix extracts after revising the window size of image subject to registration, thereby the size of adaptively correcting analysis window.

The advantage of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention is:

(1) the present invention can be implemented in the high precision autoregistration between remote sensing images under the big rotation zoom mode, has strengthened the robustness that remote sensing image registration is handled, and has improved the treatment effeciency of remotely-sensed data, for the automatic processing that realizes remotely-sensed data lays the foundation.

(2) the present invention combines angle point feature, local segmentation image and local edge contour, realizes the automatic extraction of same place jointly as local feature;

(3) the present invention introduces registration process with local configuration covariance matrix, realize the automatic extraction of rotation zoom factor between reference picture and image subject to registration, simultaneously reduce to rotate the influence of zoom factor, improve the robustness of registration process registration process in conjunction with local window self-adaptation selecting technology.

Description of drawings

Fig. 1 is based on the process flow diagram of the registration process method of characteristics of image;

Fig. 2 is the flow chart of steps of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention;

Fig. 3 is the angle point characteristic extraction step process flow diagram of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention;

Fig. 4 is the Edge Gradient Feature flow chart of steps of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention;

Fig. 5 is the optimum window size extraction step process flow diagram of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention;

Fig. 6 is the angle point characteristic matching treatment step process flow diagram of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention;

Fig. 7 a is the rectangular coordinate system synoptic diagram of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention;

Fig. 7 b is that the intrinsic coordinates of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention is a synoptic diagram;

Fig. 8 a is the change curve of coefficient of rotary at the different angle points place of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention;

Fig. 8 b is the change curve of zoom factor at the different angle points place of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention;

Fig. 9 is the analysis result figure that the optimum window size of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention extracts;

Figure 10 a is the change curve of the intersegmental similarity Dkl in the partial correction front of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention;

Figure 10 b is the change curve of the intersegmental similarity Dkl in the partial correction front of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention;

Figure 11 is the registration process figure as a result of a kind of remote sensing image registration method based on local configuration covariance matrix of the present invention.

Embodiment

The present invention is described in further detail below in conjunction with drawings and Examples.

The present invention is primarily aimed at and is difficult to set up this difficult point of mapping relations under the big rotation zoom mode between the remote sensing images feature, a kind of remote sensing image registration disposal route based on local configuration covariance matrix has been proposed, this disposal route can realize the automatic extraction of rotation zoom factor between remote sensing images, increased the matching precision between the remote sensing images feature effectively, strengthen the robustness that remote sensing image registration is handled, realized the accurate registration between remote sensing images under the big rotation zoom mode under the no artificial participation condition.

This method is at first extracted the angle point feature in the input picture, be characterized as reference with angle point, extract near the local configuration feature of angle point, and the covariance matrix that utilizes the local configuration feature extracts the rotation zoom factor between input picture, reduce to rotate the influence of zoom factor, improve the extraction precision of same place characteristic matching; Then, utilize detected same place to carry out the affine model parameter estimation; At last, carry out image resampling according to the transformation model parameter that estimates, the remote sensing images behind the acquisition registration are right.On the whole, this treatment scheme automatic extraction, optimum window size of mainly comprising the selection of feature and extraction, rotation zoom factor choose, angle point characteristic matching, affine model parameter estimation and six parts of image resampling.

A kind of remote sensing image registration method based on local configuration covariance matrix of this method as shown in Figure 2, comprises the steps:

Step 1: the selection of feature and extraction.

The selection of feature is the prerequisite that remote sensing image registration is handled with extracting, and it directly has influence on the effect of registration process, even causes realizing matching treatment.Say in principle, complete edge contour is a kind of good characteristics of image, and it not only has good station-keeping ability, and can set up mapping relations easily, therefore, for optical imagery, utilize edge contour to can be good at realizing the image registration processing.Yet for radar image, strong speckle noise makes and is difficult to extract complete edge contour, at this moment, can't obtain result preferably if directly adopt edge feature to carry out image registration.In order effectively to address this problem, adopt in the present invention angle point feature, local segmentation image are combined with the local edge contour feature, make full use of the feature that the good and local configuration feature of angle point feature location characteristic is easy to set up mapping relations, consider that simultaneously the angle point feature is usually located at the place of gray scale sudden change, therefore can effectively reduce local configuration Feature Extraction difficulty, improve the robustness of registration process.

Angle point Feature Extraction: adopt the Harris Corner Detection Algorithm to realize the image Corner Detection in this article.This algorithm judges whether to exist angle point by the gradient of gradation of image on all directions in the statistics mobile search window.If the grey scale change that gradation of image is gone up in any direction in the mobile search window is all very little, then do not comprise the angle point feature in the region of search, if the tangible grey scale change of all appearance that gradation of image is gone up in any direction in the mobile search window, then the feature in the region of search is an angle point.For effectively judging the variation of gradation of image, define angular response function C RF, it comprises two kinds of patterns and is respectively Harris pattern and Nobel pattern:

$\mathrm{CRF}=\left\{\begin{array}{cc}\det M-k\·\mathrm{trace}\left(M\right)& \left(\mathrm{Harris}\right)\\ \frac{\mathrm{trace}\left(M\right)}{\det M}& \left(\mathrm{Nobel}\right)\end{array}\right.---\left(1\right)$

Wherein: M represents gradient matrix, and its expression formula is:

$M=G\⊗\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]---\left(2\right)$

I _xGradient on the expression x direction; I _yGradient on the expression y direction; G is Gauss's template;

The expression convolution algorithm; The determinant of det representing matrix; The mark of trace representing matrix; K is a constant coefficient, generally gets 0.04.

Whether the CRF value of judging every bit is the interior maximal value (Nobel angular response function correspondence minimum value) of local neighborhood.If then this point is set to alternative angle point; If not, show that then this point is not an angle point.Utilize threshold processing to extract final angle point at last.

Fig. 3 has provided angle point and has extracted process flow diagram.

During concrete operations, use masterplate respectively $\left[\begin{array}{ccc}-1& 0& 1\\ -1& 0& 1\\ -1& 0& 1\end{array}\right]$ And masterplate $\left[\begin{array}{ccc}-1& -1& -1\\ 0& 0& 0\\ 1& 1& 1\end{array}\right]$ In orientation, image upper edge to distance to moving, calculate each locational Grad, obtain the gradient image on the both direction, and calculate the gradient product of each position correspondence, obtained three width of cloth figure so in fact, its pixel value is corresponding respectively I _x, I _yAnd I _xI _yUtilize Gauss's masterplate to carry out process of convolution with this three width of cloth image respectively then, then, calculate angular response function (CRF) value of each point respectively according to formula (1), and whether the CRF value of judging every bit is the maximal value of (as 5 * 5) in its local neighborhood, if not show this point is not angle point, if then this point is set to alternative angle point, utilize threshold processing to extract the optimum point of some at last as last result.

Step 2: local configuration extracts: on the basis of finishing the angle point feature extraction, further extract the local configuration feature of each angle point.According to the definition of angle point and angle point Feature Extraction method as can be known, the angle point feature is usually located at the place of gray scale sudden change, therefore, regional area at angle point exists tangible grey scale change, therefore, utilize simple image partition method just topography can be carried out and cause division, and then extract local configuration.The inter-class variance automatic threshold method of Cai Yonging is carried out Threshold Segmentation in the present invention.

As shown in Figure 4, according to the tonal range of image 0,1 ..., l-1} selects threshold value t that it is divided into two classes: C successively ₀: 0,1 ..., t} and C ₁: t+1, t+2 ..., l-1}.At first, the histogram of statistics topography is to image histogram normalization, and the probability distribution of each gray-scale value is:

p _i＝n _i/N，p _i≥0， $\underset{i=0}{\overset{l-1}{\mathrm{\Σ}}}{p}_i=1---\left(3\right)$

Wherein: n _iThe pixel number of representing i gray-scale value, N represents total pixel number, p _iThe shared probability of pixel of representing this gray-scale value.

At this moment, can obtain the inter-class variance σ of two classes _B ²For:

${\mathrm{\&sigma;}}_{\mathrm{<figure-callout\; id="2"\; label="B"\; filenames="S2008101023293E00072.png,S2008101023293E00081.png"\; state="\{\{state\}\}">B</figure-callout>}}^2={\mathrm{\&omega;}}_0{({\mathrm{\&mu;}}_0-{\mathrm{\&mu;}}_T)}^2+{\mathrm{\&omega;}}_1{({\mathrm{\&mu;}}_1-{\mathrm{\&mu;}}_T)}^2={\mathrm{\&omega;}}_0{\mathrm{\&omega;}}_1{({\mathrm{\&mu;}}_1-{\mathrm{\&mu;}}_0)}^2---\left(4\right)$

Wherein: ω ₀And ω ₁Represent C respectively ₀Class and C ₁The probability of occurrence of class, its computing formula is:

$\left\{\begin{array}{c}{\mathrm{\&omega;}}_0=\underset{i=0}{\overset{t}{\mathrm{\&Sigma;}}}{p}_i\\ {\mathrm{\&omega;}}_1=\underset{i=t+1}{\overset{l-1}{\mathrm{\&Sigma;}}}{p}_i\end{array}\right.---\left(5\right)$

μ ₀, μ ₁And μ _TRepresent C respectively ₀Class, C ₁The average on class and whole gray scales rank, its computing formula is:

$\left\{\begin{array}{c}{\mathrm{\&mu;}}_0=\underset{i=0}{\overset{t}{\mathrm{\&Sigma;}}}\frac{{\mathrm{ip}}_i}{{\mathrm{\&omega;}}_0}\\ {\mathrm{\&mu;}}_1=\underset{i=t+1}{\overset{l-1}{\mathrm{\&Sigma;}}}\frac{{\mathrm{ip}}_i}{{\mathrm{\&omega;}}_1}\\ {\mathrm{\&mu;}}_T=\underset{i=0}{\overset{l-1}{\mathrm{\&Sigma;}}}{\mathrm{ip}}_i\end{array}\right.---\left(6\right)$

Set different segmentation thresholds, calculate the inter-class variance of two class images under various Threshold Segmentation situations respectively, and be criterion to the maximum, can obtain optimal threshold t with inter-class variance according to formula (4), formula (5) and formula (6) ^*For:

$t^{*}\text{=}\underset{0\&le;t\&le;l-1}{\mathrm{<figure-callout\; id="0"\; label="Max"\; filenames="S2008101023293E00082.png,S2008101023293E00101.png"\; state="\{\{state\}\}">Max</figure-callout>}}{\mathrm{\&sigma;}}_{B---\left(7\right)}^2$

And utilize this threshold value to finish the two-value of image is cut apart, pixel value by each point in the neighborhood around each point in the comparative analysis bianry image and its changes and extracts the edge feature of protecting local neighborhood at last, when existing pixel value to change, show and have edge feature, when not existing pixel value to change, show and do not have edge feature, at last, extract the nearest edge feature of digression point feature as the local configuration feature by the relative position relation of analyzing angle point and edge.

Step 3: the automatic extraction of rotation zoom factor

When having bigger rotation scaled error between input picture, the local similar between two images sharply reduces with the increase of rotation scaled error, and this makes characteristic matching become difficult more.Rotate zoom factor between input picture and utilize local covariance matrix between input picture to extract effectively, and then reduce to rotate the influence of zoom factor, improve the robustness of registration process characteristic matching.

Any piece image himself all exists such intrinsic coordinates, and promptly by this coordinate system, the covariance matrix of this image is a diagonal matrix, and the eigenwert of covariance matrix has then reflected the energy distribution in the image.With the star target is that example illustrates intrinsic coordinates system, shown in Fig. 7 a, in rectangular coordinate system, there is a star target, at this moment, intrinsic coordinates for this star target is shown in Fig. 7 b, image is about the axis symmetry in intrinsic coordinates system, and its covariance matrix is a diagonal matrix, and eigenwert has reflected the distribution of image energy.The relation of the intrinsic coordinates system by contrast input picture angle point place local configuration can be extracted the twiddle factor between registering images and image subject to registration, and relatively the eigenwert of covariance matrix can obtain the zoom factor between two images.Below derive the in detail extracting method of rotation zoom factor.

Suppose I _iAnd T _jBe respectively the corresponding angle point in reference picture and the image subject to registration, then reference picture angle point I _iLocal neighborhood I ₀With image angle point T subject to registration _jLocal neighborhood T ₀Between have affined transformation f, make between two images and to satisfy:

T ₀＝f(I ₀)＝AI ₀(8)

Wherein: A is 2 * 2 affine transformation matrix.

Distributed according to local gray level earlier before extracting local yardstick twiddle factor and take Threshold Segmentation to handle to topography, gray-scale value is set to 1 greater than the point of threshold value, is set to 0 less than the point of threshold value, forms characteristic image, and note is made I and T.Though consider that the gray difference of multi-mode remote sensing images is very big, the basic structural feature of each image is consistent, therefore, still satisfy affine transformation relationship between its characteristic image.

Employing is based on the rotation zoom factor extraction method of local configuration covariance matrix.The bianry image that utilizes Threshold Segmentation to obtain calculates the covariance matrix ∑ of reference picture I and image T subject to registration respectively _IAnd ∑ _T, satisfy between two covariance matrixes: the covariance matrix ∑ that calculates reference picture I and image T subject to registration respectively _IAnd ∑ _T, satisfy between two covariance matrixes:

∑ _T＝A∑ _IA′(9)

Wherein: ${\mathrm{\&Sigma;}}_I=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{r}}_{\mathrm{ij}}{\stackrel{\&RightArrow;}{r}}_{\mathrm{ij}}^{\&prime;};$ ${\mathrm{\&Sigma;}}_T=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{r}}_{\mathrm{tj}}{\stackrel{\&RightArrow;}{r}}_{\mathrm{tj}}^{\&prime;};$ N is the pixel number of O; ${\stackrel{\&RightArrow;}{r}}_j=\left[\begin{array}{c}x\\ y\end{array}\right]$ Be the position vector of each point, subscript ' expression transposition.

Because oblique variance matrix ∑ _IAnd ∑ _TBe normal matrix, have orthogonal matrix U ₁, U ₂With diagonal matrix Λ ₁, Λ ₂Make:

${\mathrm{\&Sigma;}}_I=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{r}}_{\mathrm{ij}}{\stackrel{\&RightArrow;}{r}}_{\mathrm{ij}}^{\&prime;}=U_1{\mathrm{\&Lambda;}}_1{U}_1^{\&prime;}$ ${\mathrm{\&Sigma;}}_T=\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\stackrel{\&RightArrow;}{r}}_{\mathrm{tj}}{\stackrel{\&RightArrow;}{r}}_{\mathrm{tj}}^{\&prime;}={\mathrm{<figure-callout\; id="2"\; label="U"\; filenames="S2008101023293E00072.png,S2008101023293E00081.png"\; state="\{\{state\}\}">U</figure-callout>}}_2{\mathrm{\&Lambda;}}_2{\mathrm{<figure-callout\; id="2"\; label="U"\; filenames="S2008101023293E00072.png,S2008101023293E00081.png"\; state="\{\{state\}\}">U</figure-callout>}}_2^{\&prime;}---\left(10\right)$

Wherein: $U_1=\left[\begin{array}{cc}\cos {\mathrm{\&theta;}}_1& \sin {\mathrm{\&theta;}}_1\\ -\sin {\mathrm{\&theta;}}_1& \cos {\mathrm{\&theta;}}_1\end{array}\right];$ ${\mathrm{<figure-callout\; id="2"\; label="U"\; filenames="S2008101023293E00072.png,S2008101023293E00081.png"\; state="\{\{state\}\}">U</figure-callout>}}_2=\left[\begin{array}{cc}\cos {\mathrm{\&theta;}}_2& \sin {\mathrm{\&theta;}}_2\\ -\sin {\mathrm{\&theta;}}_2& \cos {\mathrm{\&theta;}}_2\end{array}\right];$

${\mathrm{\&Lambda;}}_1=\left[\begin{array}{cc}{\mathrm{\&lambda;}}_{11}& 0\\ 0& {\mathrm{\&lambda;}}_{12}\end{array}\right],$ λ ₁₁＞λ ₁₂； ${\mathrm{\&Lambda;}}_2=\left[\begin{array}{cc}{\mathrm{\&lambda;}}_{21}& 0\\ 0& {\mathrm{\&lambda;}}_{22}\end{array}\right],$ λ ₂₁＞λ ₂₂。

Following formula substitution formula (9) is had:

U ₂Λ ₂U′ ₂＝AU ₁Λ ₁U′ ₁A′(11)

Order $B_{}$ ${}_2={\mathrm{<figure-callout\; id="2"\; label="U"\; filenames="S2008101023293E00072.png,S2008101023293E00081.png"\; state="\{\{state\}\}">U</figure-callout>}}_2{\mathrm{\&Lambda;}}_2^{0.5},$ $B_1=U_1{\mathrm{\&Lambda;}}_1^{0.5}$ Have:

$A=B_2{B}_1^{-1}---\left(12\right)$

Like this, distribute according to the partial structurtes of reference picture and image subject to registration and can extract twiddle factor between image

With zoom factor λ.

λ _i＝(λ _2i/λ _1i)(13)

When specific implementation rotation zoom factor extracts, need the corresponding relation between earlier definite characteristic image, when determining the corresponding relation between image, dark, bright image in bright, dark image in the corresponding image subject to registration of bright in the comparative analysis reference picture, dark image and the reference picture in the corresponding image subject to registration of bright, dark image, the ratio of anglec of rotation difference and zoom factor under the both of these case, it is littler to choose anglec of rotation difference, diaxon rotation zoom factor λ ₁And λ ₂The more approaching corresponding relation of ratio as the corresponding relation between two images.

Shown in Fig. 8 a, be the extraction result of the anglec of rotation of the different angle points place of input picture in this example estimation, the zoom factor that Fig. 8 b has provided each point extracts the result.From scheming as seen, the variation range of the twiddle factor of extraction is in 1.2 degree, and the variation range of zoom factor is in 0.14, and they all have very high consistance separately, has reacted the rotation convergent-divergent relation between input picture.

Step 4: optimum window size is chosen

When having bigger rotation and scale error between two width of cloth images, the test performance of local matching criterior commonly used such as related coefficient, crosscorrelation coefficient, interactive information etc. sharply descends, and the reason that causes test performance to descend mainly contains two aspects:

1: because the test performance that twiddle factor and zoom factor cause descends;

2, owing to the influence of zoom factor, under identical window size, the difference of picture material causes test performance to descend.Therefore, when having bigger rotation and scaled error between two width of cloth images, need extract twiddle factor and zoom factor between two width of cloth images according to the content of image, the size of adaptive adjustment analysis window is eliminated the influence of rotation scaled error.

The cardinal rule that self-adaptation is adjusted the analysis window size is as follows: when selected window size is optimum window size, picture material unanimity in two windows, the ratio of two image local configuration covariance matrix eigenwerts has reflected the zoom factor between two images at this moment, the ratio of analysis window has reflected the zoom factor between two images equally in while two images, in theory will, when reaching optimum window size, two ratios are consistent, therefore, can choose the analysis window size to the zoom factor λ that local configuration covariance matrix extracts by the ratio p and the input picture of selected window size in comparative analysis reference picture and the image subject to registration.As shown in Figure 5, specific analytical method is as follows:

The zoom factor λ that the statistics local configuration covariance matrix extracts, calculate the ratio p that twiddle factor and zoom factor calculate selected window size in comparative analysis reference picture and the image subject to registration then, and choose the difference of λ and p under the different windows size condition in the comparative analysis image subject to registration, if it is minimum that the difference of λ and p reaches, then obtain the size of the size of this window subject to registration as optimum window; Otherwise recomputate the zoom factor λ that local configuration covariance matrix extracts after revising the window size of image subject to registration, thereby the size of adaptively correcting analysis window.

If the interior partial image texture characteristic of local window is identical in the topography in the reference picture in the local window and the image subject to registration, at this moment, the ratio of selected window size has also just in time reflected the zoom factor between two images in reference picture and the image subject to registration, therefore, will be consistent according to the zoom factor that ratio calculated of selected window size in reference picture and the image subject to registration with the zoom factor that is extracted according to formula (13).

If the partial image texture characteristic in topography in the reference picture in the local window and the image subject to registration in the local window is inequality, the zoom factor that is extracted according to the zoom factor that ratio obtained of selected window size in reference picture and the image subject to registration with according to formula (13) all can not be correct extraction topography between zoom factor, at this moment, will there be certain deviation between two results.

Therefore, can judge whether to obtain optimum window size, promptly hour reach optimum window size when difference reaches by the difference of the zoom factor that extracted under two kinds of situations relatively.Fig. 9 has provided the analysis result that optimum window size extracts, there is the optimum window of selecting between 23～25, to fluctuate as seen from the figure, by final registration process as can be known the theoretical optimum window between two images be 23, have the consistance of height by the relatively more visible optimum window size that extracts, well reacted the rotation convergent-divergent relation between input picture.

Step 5: angle point characteristic matching

Utilize the local edge contour feature to realize the coupling of angle point feature.At first utilize the rotation zoom factor and the optimum window size that extract that the edge contour that extracts in the step 1 is resampled, reduce to rotate the influence of zoom factor characteristic matching.

Then, according to the spatial relation between each point the local edge feature after resampling is carried out the Freeman coding, the similarity measurement formula of encoding between the line segment that back two length are N through curve is shown below:

$D_{\mathrm{kl}}=\frac{1}{N}\underset{j=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\cos \frac{\mathrm{\&pi;}}{4}(a_{k+j}-b_{l+j})---\left(14\right)$

Wherein: a _K+jAnd b _L+jRepresent that respectively two line segments carry out the resulting code value in Freeman coding back; D _KlBe the value between 0 to 1, two line segments are coupling more, and it is big more that its value levels off to 1, two straight line difference more, and its value levels off to 0, two line segment when mating fully approximately, and its value is 1.

Add up the similarity of two line segments, the similarity D when between two line segments _KlWhen satisfying following condition, think that this is a same place to angle point.

$D_{\mathrm{kl}}=\underset{k\&Element;A}{\max }\left(D_{\mathrm{kl}}\right)$ $D_{\mathrm{kl}}=\underset{l\&Element;B}{\max }\left(D_{\mathrm{kl}}\right)$ D _kl＞D _th(15)

Wherein: A represents the angle point collection that extracts in the reference picture; B represents the angle point collection that extracts in the image subject to registration; D _ThThe threshold value of expression curve similarity.

So far, can realize the thick coupling of angle point feature between remote sensing images by above-mentioned steps, mapping relations have between points been set up, yet because the difference between noise and gradation of image, inevitably exist error matching points right, in order to improve the precision of registration process, be necessary that the match point of further rejecting mistake is right on the basis of thick coupling.

For this reason, can further utilize the mutual support characteristic between relative position relation between point set and point set to reject error matching points.At first, top described initial matching point is sorted to the size according to the similarity test.Consider the space structure that 4 pairs of reference mark can fine reaction reference mark part, and the reference mark group that is made of 4 pairs of reference mark, the reference mark group that sequential search is all are carried out being combined into according to putting in order in all reference mark.Suppose to satisfy between them the affined transformation model, when these 4 pairs of reference mark when spatially not having the conllinear phenomenon, very little according to these 4 pairs of reference mark to the affine matrix and the error between actual affine matrix of estimation, affined transformation is carried out at the reference mark that the affine matrix of utilization estimation is treated in the registering images, and with 4 pairs of reference mark to average position error be that index weighs whether be that correct reference mark is right, when the reference mark to average position error during less than preset threshold, can think that to find correct reference mark right, stop search.When the average position error at reference mark greater than preset threshold or the 4 pairs of reference mark when spatially having the conllinear phenomenon, then think the mistake match point, continue search, all analyze up to the reference mark of initial extraction and finish.

Determining on the right basis, the 4 pairs of correct reference mark, testing remaining reference mark to whether being that correct reference mark is right.Order is right to introducing correct reference mark with reference mark to be tested, according to correct reference mark to estimating affine matrix, recomputate the reference mark to average position error, when average position error during greater than setting threshold, think that initiate reference mark is right to being wrong reference mark, it is deleted from initial matching point centering.Otherwise, think the reference mark of new introducing to right for correct reference mark, it is right to continue all the other reference mark of search, and it is right to extract all correct reference mark.Fig. 6 provides the process flow diagram that the angle point characteristic matching is handled.

Shown in Figure 10 a, 10b, Figure 10 a is the similarity D that does not extract under rotation zoom factor and the optimum window size _KlChange curve, Figure 10 b is for extracting the similarity D behind rotation zoom factor and the optimum window size _KlChange curve.From the result, as seen,, improved the robustness of angle point characteristic matching greatly owing to can eliminate the influence of rotation zoom factor between topography.

Step 6: model parameter estimation and image resampling

After finishing Feature Points Matching, can utilize the parameter in the least square method estimation affine model.Concrete method of estimation is as follows:

Suppose to have extracted altogether between input picture the N in the step 4 to the reference mark, be respectively { ((x _R1, y _R1), (x _W1, y _W1)), ((x _R2, y _R2), (x _W2, y _W2)) ..., ((x _RN, y _RN), (x _WN, y _WN)), consider and satisfy affine relation between image that at this moment, arbitrarily the reference mark is to satisfying following relational expression:

$\left[\begin{array}{c}{x}_{\mathrm{wi}}\\ y_{\mathrm{wi}}\end{array}\right]=\left[\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{ri}}\\ y_{\mathrm{ri}}\end{array}\right]+\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;y}\end{array}\right]=\left[\begin{array}{cccccc}1& 0& x_{\mathrm{<figure-callout\; id="0"\; label="ri\; y\; ri"\; filenames="S2008101023293E00082.png,S2008101023293E00101.png"\; state="\{\{state\}\}">ri</figure-callout>}}& y_{\mathrm{ri}}{}_{}& 0& 0\\ 0& 1& 0& 0& x_{\mathrm{ri}}& y_{\mathrm{ri}}\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;y}\\ a_{11}\\ a_{12}\\ a_{21}\\ a_{22}\end{array}\right]---\left(16\right)$

To extract the equal substitution equation in all reference mark, can obtain following system of equations:

$Y=\left[\begin{array}{c}{x}_{w1}\\ y_{\mathrm{wl}}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ x_{\mathrm{wN}}\\ y_{\mathrm{wN}}\end{array}\right]=\left[\begin{array}{cccccc}1& 0& x_{r1}& y_{r1}& 0& 0\\ 0& 1& 0& 0& x_{r1}& y_{r1}\\ & & & \&CenterDot;& & \\ & & & \&CenterDot;& & \\ & & & \&CenterDot;& & \\ 1& 0& x_{\mathrm{<figure-callout\; id="0"\; label="rN\; y\; rN"\; filenames="S2008101023293E00082.png,S2008101023293E00101.png"\; state="\{\{state\}\}">rN</figure-callout>}}& y_{\mathrm{rN}}{}_{}& 0& 0\\ 0& 1& 0& 0& x_{\mathrm{rN}}& y_{\mathrm{rN}}\end{array}\right]=\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;y}\\ a_{11}\\ a_{12}\\ a_{21}\\ a_{22}\end{array}\right]=H\&CenterDot;\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;y}\\ a_{11}\\ a_{12}\\ a_{21}\\ a_{22}\end{array}\right]---\left(17\right)$

After each position coordinates substitution system of equations to the reference mark, the estimated value that can obtain affine parameter is shown below:

$\left[\begin{array}{c}\mathrm{\&Delta;x}\\ \mathrm{\&Delta;y}\\ a_{11}\\ a_{12}\\ a_{21}\\ a_{22}\end{array}\right]={(H^T\&CenterDot;H)}^{-1}\&CenterDot;H^T\&CenterDot;Y---\left(18\right)$

Wherein: subscript T represents transposition, and subscript-1 expression is inverted.

After determining transformation parameter, can be according to affine parameter to input picture to processings that resample, the image that can obtain coupling mutually is right, adopts the resampling of bilinear interpolation method realization image in this Processing Algorithm.

At first, calculate matrix of coefficients H and matrix H in the formula (17) according to the position of each point in the reference diagram ^TH is then with matrix H ^T.H be decomposed into nonsingular diagonal matrix D and nonsingular triangle row matrix Λ down,

$D=\left[\begin{array}{cccccc}{d}_1& 0& 0& 0& 0& 0\\ 0& {\mathrm{<figure-callout\; id="2"\; label="d"\; filenames="S2008101023293E00072.png,S2008101023293E00081.png"\; state="\{\{state\}\}">d</figure-callout>}}_2& 0& 0& 0& 0\\ 0& 0& d_3& 0& 0& 0\\ 0& 0& 0& {\mathrm{<figure-callout\; id="4"\; label="d"\; filenames="S2008101023293E00082.png,S2008101023293E00091.png"\; state="\{\{state\}\}">d</figure-callout>}}_4& 0& 0\\ 0& 0& 0& 0& d_5& 0\\ 0& 0& 0& 0& 0& {\mathrm{<figure-callout\; id="6"\; label="d"\; filenames="S2008101023293E00072.png,S2008101023293E00081.png"\; state="\{\{state\}\}">d</figure-callout>}}_6\end{array}\right],$ $\mathrm{\&Lambda;}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ {\mathrm{\&lambda;}}_{21}& 1& 0& 0& 0& 0\\ {\mathrm{\&lambda;}}_{31}& {\mathrm{\&lambda;}}_{32}& 1& 0& 0& 0\\ {\mathrm{\&lambda;}}_{41}& {\mathrm{\&lambda;}}_{42}& {\mathrm{\&lambda;}}_{43}& 1& 0& 0\\ {\mathrm{\&lambda;}}_{51}& {\mathrm{\&lambda;}}_{52}& {\mathrm{\&lambda;}}_{53}& {\mathrm{\&lambda;}}_{54}& 1& 0\\ {\mathrm{\&lambda;}}_{61}& {\mathrm{\&lambda;}}_{62}& {\mathrm{\&lambda;}}_{63}& {\mathrm{\&lambda;}}_{64}& {\mathrm{\&lambda;}}_{65}& 1\end{array}\right]---\left(19\right)$

Wherein: $x_{\mathrm{ij}}=h_{\mathrm{ij}}-\underset{k=1}{\overset{j-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{ik}}{\mathrm{\&lambda;}}_{\mathrm{jk}},$ ${\mathrm{\&lambda;}}_{\mathrm{ij}}=\frac{x_{\mathrm{ij}}}{d_j},$ $d_i=h_{\mathrm{ii}}-\underset{k=1}{\overset{j-1}{\mathrm{\&Sigma;}}}{x}_{\mathrm{ik}}{\mathrm{\&lambda;}}_{\mathrm{jk}}\&CenterDot;$

Ask for the inverse matrix D of matrix D and matrix Λ respectively ^-1And Λ ^-1

$D^{-1}=\left[\begin{array}{cccccc}\frac{1}{d_1}& 0& 0& 0& 0& 0\\ 0& \frac{1}{{\mathrm{<figure-callout\; id="2"\; label="d"\; filenames="S2008101023293E00072.png,S2008101023293E00081.png"\; state="\{\{state\}\}">d</figure-callout>}}_2}& 0& 0& 0& 0\\ 0& 0& \frac{1}{d_3}& 0& 0& 0\\ 0& 0& 0& \frac{1}{{\mathrm{<figure-callout\; id="4"\; label="d"\; filenames="S2008101023293E00082.png,S2008101023293E00091.png"\; state="\{\{state\}\}">d</figure-callout>}}_4}& 0& 0\\ 0& 0& 0& 0& \frac{1}{d_5}& 0\\ 0& 0& 0& 0& 0& \frac{1}{{\mathrm{<figure-callout\; id="6"\; label="d"\; filenames="S2008101023293E00072.png,S2008101023293E00081.png"\; state="\{\{state\}\}">d</figure-callout>}}_6}\end{array}\right],$ $\mathrm{\&Lambda;}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ {\mathrm{\&rho;}}_{21}& 1& 0& 0& 0& 0\\ {\mathrm{\&rho;}}_{31}& {\mathrm{\&rho;}}_{32}& 1& 0& 0& 0\\ {\mathrm{\&rho;}}_{41}& {\mathrm{\&rho;}}_{42}& {\mathrm{\&rho;}}_{43}& 1& 0& 0\\ {\mathrm{\&rho;}}_{51}& {\mathrm{\&rho;}}_{52}& {\mathrm{\&rho;}}_{53}& {\mathrm{\&rho;}}_{54}& 1& 0\\ {\mathrm{\&rho;}}_{61}& {\mathrm{\&rho;}}_{62}& {\mathrm{\&rho;}}_{63}& {\mathrm{\&rho;}}_{64}& {\mathrm{\&rho;}}_{65}& 1\end{array}\right]---\left(20\right)$

Wherein: $\left\{\begin{array}{c}{\mathrm{\&rho;}}_{\mathrm{ij}}=-\underset{k=j}{\overset{i-1}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_{\mathrm{ik}}{\mathrm{\&rho;}}_{\mathrm{kj}}\\ {\mathrm{\&rho;}}_{\mathrm{ii}}=1\end{array}\right.(j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n-1;i=j+1,j+2,\&CenterDot;\&CenterDot;\&CenterDot;,n)$

Like this, utilize inverse matrix D ^-1And Λ ^-1Can ask for matrix H ^TThe inverse matrix of H, right position, reference mark can calculate affine transformation parameter in associate(d) matrix H and the image subject to registration.

Table 1 has provided the registration process result, Figure 11 has provided the change curve of the root-mean-square error of position, reference mark, the position at reference mark in reference picture X coordinate and the Y coordinate representation reference picture wherein, the position at reference mark in image X coordinate subject to registration and the Y coordinate representation image subject to registration, the room parameter that back X coordinate and Y coordinate representation utilization extract of resampling is treated reference mark in the registering images and is carried out the position coordinates that obtains after the conversion process, the distance that departs between the position, reference mark in the position at reference mark and the reference picture after the root-mean-square error of position is represented to resample, from table 1 and Figure 11 as can be seen, there are 25 degree rotation errors, 1.4 under the zoom factor condition doubly, the 10 pairs of same places are 0.94 pixel cell to the maximum through the root mean square of position after the treatment for correcting, can realize that registration error is better than 1 pixel cell.

Table 1 registration process result

Sequence number	Reference picture X coordinate	Reference picture Y coordinate	Image X coordinate subject to registration	Image Y coordinate subject to registration	Back X coordinate resamples	Back Y coordinate resamples	The root-mean-square error of position
								0	652	856	943	564	651.78	855.95	0.23
1	620	679	989	390	619.96	678.87	0.14
								2	498	655	889	317	498.54	654.95	0.54
3	188	520	665	64	188.78	520.30	0.84
								4	298	758	664	325	298.11	757.20	0.81
5	45	874	386	324	45.88	873.67	0.94

6	316	782	670	355	316.22	781.85	0.27
								7	104	549	576	54	103.94	548.82	0.19
8	442	442	928	100	442.21	441.87	0.25
								9	447	591	869	237	446.64	590.91	0.37

Claims (4)

1. the remote sensing image registration method based on local configuration covariance matrix is characterized in that, comprises the steps:

Step 1: from each image channel, read in each image simultaneously, adopt the Harris Corner Detection Algorithm to realize the image Corner Detection, judge whether to exist angle point by the angular response functional value of adding up each point in the mobile search window;

If the angular response function C RF value of any is the maximal value in the local neighborhood, then this point is set to alternative angle point; Otherwise show that this point is not an angle point; After in having judged image, having a few, utilize threshold processing to extract final angle point;

Step 2: adopt inter-class variance automatic threshold method that reference picture and image subject to registration are carried out local threshold and cut apart;

Utilize the optimal threshold that extracts that topography is cut apart, obtain local bianry image, and obtain the local edge profile according to the local bianry image that extracts;

Step 3: the automatic extraction of rotation zoom factor; The automatic extraction of described rotation zoom factor is based on the local configuration covariance matrix realization, concrete steps comprise: the bianry image that utilizes Threshold Segmentation to obtain, calculate the covariance matrix of reference picture and image subject to registration respectively, distribute according to the partial structurtes of reference picture and image subject to registration and can extract twiddle factor and zoom factor between image;

Step 4: choose the analysis window size by the ratio of selected window size in comparative analysis reference picture and the image subject to registration and the zoom factor that utilizes input picture to extract, extract optimum window;

Step 5: utilize the rotation zoom factor and the optimum window size that extract that the edge contour that extracts in the step 2 is resampled, reduce to rotate the influence of zoom factor to characteristic matching; And the edge feature after resampling carried out the Freeman coding, and add up similarity between two line segments by the chain code correlation technique; Utilize the mutual support characteristic between relative position relation between point set and point set to reject error matching points;

Step 6: after finishing Feature Points Matching, utilize the parameter in the least square method estimation affine model; After determining transformation parameter, according to affine parameter to input picture to processings that resample, it is right to obtain the image that mates mutually.

2. according to the described a kind of remote sensing image registration method of claim 1, it is characterized in that based on local configuration covariance matrix:

In the described step 3, at the automatic extraction of local segmentation image stretching rotation zoom factor; It is bright, dark two class images that passing threshold is cut apart image segmentation, adds up two-part covariance matrix and rotation zoom factor then respectively, determines the corresponding relation between image, and extracts corresponding rotation zoom factor.

3. according to the described a kind of remote sensing image registration method of claim 2, it is characterized in that based on local configuration covariance matrix:

In the described step 3, when determining the corresponding relation between image, these two kinds of corresponding relations of dark, bright image in bright, dark image in the corresponding image subject to registration of bright in the comparative analysis reference picture, dark image and the reference picture in the corresponding image subject to registration of bright, dark image, choose the ratio of twiddle factor and zoom factor, it is littler to choose twiddle factor, and the more approaching corresponding relation of the ratio of diaxon rotation zoom factor is as the corresponding relation between two images.

4. according to the described a kind of remote sensing image registration method of claim 1, it is characterized in that based on local configuration covariance matrix:

In the described step 4, choose the analysis window size, extract optimum window, be specially: the zoom factor λ that the statistics local configuration covariance matrix extracts, calculate twiddle factor and zoom factor then, calculate the ratio p of selected window size in comparative analysis reference picture and the image subject to registration, and choose the difference of λ and p under the different windows size condition in the comparative analysis image subject to registration, if it is minimum that both differences reach, then obtain the size of the size of this window subject to registration as optimum window, otherwise recomputate the zoom factor that local configuration covariance matrix extracts after revising the window size of image subject to registration, thereby the size of adaptively correcting analysis window.

Images