CN1987896A - High resolution SAR image registration processing method and system - Google Patents

Abstract

This invention has provided a kind of processing method and system for the image registration of the synthetic aperture radar with high resolution. The method includes the following steps: (a) To build pyramid image and extract the mixture characteristics of the image, and then get the edge information and take it as matching input data. (b) To match the entire image rudely basing on the characteristics extracted. (c) To extract the feature corner of both the major and minor images after the rude match, and carry on the cross search for homonymic point to get the homonymic point pairs. (d) To divide the adaptive region according to the homonym point pairs. The system provided includes: the extraction cell for image data, the preprocessing cell for image and the matching cell for image. The matching precision at global sub-pixel level can be perfectly realized by this method and system. And so, it provide the base technology sustains for the further generalization of SAR data in the application fields such as map drawing, change detection, and the rebuilding of the three dimensional terrain.

Description

High resolution SAR image registration processing method and system

Technical field

The present invention relates to the image processing technique field, be specifically related to Remote Sensing Image Matching disposal route and system, especially the autoegistration method of synthetic aperture radar (SAR) image and system.

Background technology

In the remote sensing image integrated analysis, Image registration is the basis of remotely-sensed data analysis and processing, is the key link that visual fusion, dynamic change detection, three dimensional terrain reconstruction, remote sensing image such as inlay at application technology.Looked back nearly ten years, the technical development of remote sensing image autoregistration is rapid, but relevant synthetic aperture radar (SAR) image, the especially research work of high resolution SAR image autoregistration but are in the starting stage.So far, the autoregistration algorithm of more existing relevant SAR images is proposed in succession; 2000, Thepaut etc. are from the visual fusion angle, a kind of autoregistration algorithm about ERS and SPOT image has been proposed, yet this algorithm is based on the figure of the earth, earth rotation, a large amount of prioris such as how much imaging modes of satellite orbit parameter and sensor, condition is difficult to be met under many circumstances; Calendar year 2001, Dare etc. have proposed a kind of about SAR and SPOT image autoregistration model, use image edge to carry out characteristic matching, yet its processing procedure need human intervention, and are not good to the image support of big map sheet; 2003, Chen etc. proposed a kind of matching algorithm based on image MI (Mutual Information), yet this algorithm exists the processing map sheet little, the not high deficiency of registration accuracy.

Summary of the invention

Given this, for solving airborne, the automatic/semi-automatic registration problems of satellite-borne SAR image of high resolving power.The invention provides a kind of high resolution SAR image registration processing method, comprehensively adopt the method for characteristic matching and gray scale coupling, implement by thick to smart matching strategy; In the smart coupling,, not only make coupling reach the precision of overall sub-pixel-level, and increased the dirigibility and the robustness of algorithm based on the proposition of least square imagery zone self-adaptation block algorithm.This just further promotes the use of the technical support that provides the foundation for the SAR data in applications such as map making, change-detection, three dimensional terrain reconstructions.High resolution SAR image registration processing method provided by the present invention comprises step:

(a) extract the image texture feature, obtain marginal information and import data as coupling;

(b) carry out the thick coupling of the image overall situation based on the image texture feature of being extracted, obtain the main and auxiliary image behind the registration;

(c) the feature angle point of main and auxiliary image behind the thick registration of extraction, and carry out the same place intersection search, it is right to obtain same place;

(d) handle carrying out the adaptive region piecemeal according to same place.

Described same place is the unique point of corresponding same atural object on the main and auxiliary image; The same characters of ground object point of main and auxiliary image is right to constituting same place, and same place is to also can be described as reference mark (Control Points or Tie-Points) sometimes.

Wherein step (a) comprising: create pyramid image; Adopt the Canny algorithm, main and auxiliary image is carried out edge extracting; Edge image is removed noise processed; According to edge length described edge image is further handled.

Wherein adopt quick Hausdorff distance algorithm slightly to mate in the step (b).

Wherein adopt even Distribution Strategy to extract the feature angle point in the step (c);

Wherein adopting main and auxiliary image intersection search method in the step (c), is reference with main and auxiliary image feature angle point respectively, according to the transformation relation of the thick coupling of feature, seeks same place in the image field of search; Then the gained same place to giving corresponding call number; With the gained same place to merging processing.

Wherein adopt in the step (d) based on least square region adaptivity block algorithm, described adaptive region piecemeal is handled, and comprises step:

(d1) prime area piecemeal obtains initial piecemeal subimage same place to error equation (IV) according to formula (III),

x ₂＝a ₀+a ₁x+a ₂y

(III)

y ₂＝b ₀+b ₁x+b ₂y

$V_x^i=\frac{\∂x}{{\∂a}_0^i}{\mathrm{da}}_0^i+\frac{\∂x}{{\∂a}_1^i}{\mathrm{da}}_1^i+\frac{\∂x}{{\∂a}_2^i}{\mathrm{da}}_2^i-l_x^i$

$V_y^i=\frac{\∂y}{\∂b_0^i}{\mathrm{db}}_0^i+\frac{\∂y}{\∂b_1^i}{\mathrm{db}}_1^i+\frac{\∂y}{\∂b_2^i}{\mathrm{db}}_2^i-l_y^i---\left(\mathrm{IV}\right)$

$l_x^i=x_2^i-(a_0^i+a_1^{i}x+a_2^{i}y)$

$l_y^i=y_2^i-(b_0^i+b_1^{i}x+b_2^{i}y)$

(d2) iterative computation deformation parameter corrected value da ₀ ⁱ, da ₁ ⁱ, da ₂ ⁱ, db ₀ ⁱ, db ₁ ⁱ, db ₂ ⁱAnd corresponding deformation parameter a ₀ ⁱ, a ₁ ⁱ, a ₂ ⁱ, b ₀ ⁱ, b ₁ ⁱ, b ₂ ⁱ

(d3) subimage block division, rollback increases the reference mark and handles;

(d4) revise the right image coordinate value of same place;

(d5) use corresponding deformation parameter to carry out geometric transformation to the subimage block that comprises in the auxilliary image successively, obtain the image behind the registration.

According to a scheme of the present invention, a kind of high resolution SAR image registration disposal system also is provided, comprising:

The image data extraction unit is used to extract the image texture characteristic, obtains marginal information, and this information is sent to the image pretreatment unit as thick coupling input data;

The image pretreatment unit carries out the thick coupling of the overall situation to the image that extracts the image texture feature; Obtain the main and auxiliary image information behind the thick registration, and relevant information is offered the Image registration processing unit processes;

The Image registration processing unit is used to extract the feature angle point of the main and auxiliary image that obtains after the thick coupling, and carries out the same place intersection search, and it is right to obtain same place, and handles carrying out the adaptive region piecemeal according to same place.

Wherein said image data extraction unit is handled the step of image, comprising: create pyramid image; Adopt the Canny algorithm, main and auxiliary image is carried out edge extracting; Edge image is removed noise processed; According to edge length described edge image is further handled.

Wherein in the image pretreatment unit, adopt quick Hausdorff distance algorithm slightly to mate.

Wherein adopt even Distribution Strategy to extract the feature angle point in the Image registration processing unit;

Wherein adopting main and auxiliary image intersection search method in the Image registration processing unit, is reference with main and auxiliary image feature angle point respectively, according to the transformation relation of the thick coupling of feature, seeks same place in the image field of search; Then the gained same place to giving corresponding call number; With the gained same place to merging processing.

Wherein adopt in the Image registration processing unit based on least square region adaptivity block algorithm, described adaptive region piecemeal is handled, and comprises step:

(d1) prime area piecemeal obtains initial piecemeal subimage same place to error equation (IV) according to formula (III),

x ₂＝a ₀+a ₁x+a ₂y

(III)

y ₂＝b ₀+b ₁x+b ₂y

$V_x^i=\frac{\∂x}{\∂a_0^i}{\mathrm{da}}_0^i+\frac{\∂x}{\∂a_1^i}{\mathrm{da}}_1^i+\frac{\∂x}{\∂a_2^i}{\mathrm{da}}_2^i-l_x^i$

$V_y^i=\frac{\∂y}{\∂b_0^i}{\mathrm{db}}_0^i+\frac{\∂y}{\∂b_1^i}{\mathrm{db}}_1^i+\frac{\∂y}{\∂b_2^i}{\mathrm{db}}_2^i-l_y^i---\left(\mathrm{IV}\right)$

$l_x^i=x_2^i-(a_0^i+a_1^{i}x+a_2^{i}y)$

$l_y^i=y_2^i-(b_0^i+b_1^{i}x+b_2^{i}y)$

(d2) iterative computation deformation parameter corrected value da ₀ ⁱ, da ₁ ⁱ, da ₂ ⁱ, db ₀ ⁱ, db ₁ ⁱ, db ₂ ⁱAnd corresponding deformation parameter a ₀ ⁱ, a ₁ ⁱ, a ₂ ⁱ, b ₀ ⁱ, b ₁ ⁱ, b ₂ ⁱ

(d3) subimage block division, rollback increases the reference mark and handles;

(d4) revise the right image coordinate value of same place;

(d5) use corresponding deformation parameter to carry out geometric transformation to the subimage block that comprises in the auxilliary image successively, obtain the image behind the registration.

Description of drawings

Fig. 1 is the process flow diagram of high resolution SAR image registration processing method of the present invention;

Fig. 2 is a processing procedure synoptic diagram according to a particular embodiment of the invention;

The quick Hausdorff distance thick flow chart that mate of Fig. 3 for being adopted in the method for the present invention;

Fig. 4 is the flow chart of the same place intersection search that adopted in the method for the present invention;

Fig. 5 is the flow chart of the LS-MMSE algorithm flow that adopted in the method for the present invention.

Embodiment

Now the specific implementation of high resolution SAR image registration processing method of the present invention is described in detail.

Fig. 2 is a processing procedure synoptic diagram according to a particular embodiment of the invention, with reference to Fig. 2.

Texture feature extraction is for providing image texture information based on the thick coupling of Hausdorff Distance global characteristics, and as the input data of matching process.In order to obtain main marginal information, suppress noise spot, adopt following method, it comprises:

Step 1: the step of creating pyramid image;

Step 2: adopt the Canny algorithm, main and auxiliary image is carried out the step of edge extracting.It should be noted that creating the pyramid image step has reached certain image smoothing effect, therefore, in the Canny algorithm is carried out, does not need to comprise in addition Gaussian Filtering Processing process; Then use the Hysteresis method, the input edge image is carried out following processing remove noise, concrete steps are:

1. scan all marginal points, range value is greater than T ₁Obtain keeping, and be designated as correct point;

2. the marginal point range value is less than T ₂Directly remove;

3. the sweep amplitude value is between T ₁, T ₂All marginal points, if there has been correct point in eight neighborhoods of this point, then also be designated as correct point; Otherwise, directly remove.Repeat edge image is carried out this computing, end until stablizing.

Step 3: the edge image of handling via Hysteresis is carried out edge length detect step.Setting threshold T, edge length obtains keeping greater than the edge of T, otherwise, directly remove.

The quick Hausdorff distance thick flow chart that mate of Fig. 3 for being adopted in the method for the present invention is with reference to Fig. 3.

Given bidimensional characteristic sequence group A and B, the Hausdorff distance h ^f(B A) is defined as follows:

$h^f(B,A)\stackrel{\mathrm{def}}{=}\mathrm{fth}\min \left|\right|b-a\left|\right|---\left(I\right)$

In the formula, The fractile of expression g (x) in sequence X, f ∈ (0,1).For example, the 1th representative is worth most, and 1/2th represents intermediate value.‖ b-a ‖ represents an a, the Euclidean distance between the b.

What the present invention adopted is briefly described below based on the Hausdorff distance algorithm:

Use fast algorithm, calculate main image E ₁Range conversion.E ₁Range conversion is defined as follows: $\mathrm{\Δ}(x,y)=\underset{i\∈E_1}{\min }\left|\right|\left(\mathrm{xy}\right)-i\left|\right|$

The storehouse (Stack) of structure storage unit R.The purpose of the thick coupling of global characteristics is the best translation transformation t[T that seeks between main and auxiliary image _x, T _y], promptly satisfy relational expression:

h ^f(t(E ₂)，E ₁)≤τ+ε (II)

In the formula, E ₁, E ₂The pyramid edge image of expression after via feature extraction; τ is a pre-set threshold, t (E ₂) represent E ₂Carry out the t conversion, ε=.R represents the two-dimensional space of translation transformation.For given R, t ^lExpression T _x, T _yMinimalization T all _x ^l, T _y ^lThe time, i.e. t ^l(T _x ^l, T _y ^l) situation; In like manner, t ^hExpression T _x, T _yAll get maximum value T _x ^h, T _y ^hThe time, i.e. t ^h(T _x ^h, T _y ^h) situation.For marginal point P ∈ E arbitrarily ₂, t ∈ R, conversion t (P) must fall within E ₁In corresponding certain rectangular area.This zone is defined as the indefinite zone of P.Therefore, for P ∈ E ₂, the just corresponding indefinite zone of each point.Initialization storehouse (Stack) makes storehouse (Stack) comprise an initial cell R, makes that the R initial grade is a zero level.It should be noted that the corresponding indefinite zone of initial cell R should select enough greatly, make it can comprise Target Transformation function t, satisfy relational expression (II).

Seek best translation transformation amount t.

A) the translation transformation difference and the maximum indefinite zone of calculating unit R correspondence to be selected

Translation transformation difference: Δ t=t ^h-t ^l

Indefinite zone: $\underset{(x_i,y_i)\∈E_2}{\max }(T_{x_i},T_{y_i})$

B) the division unit R is quartern subelement R _i(i=1,2,3,4), R _iCorresponding grade adds one for his father's cell level number, R _iBe designated as respectively: { (T _x ^l, T _y ^l), (T _x ^m, T _y ^m), { (T _x ^m, T _y ^l), (T _x ^h, T _y ^m), { (T _x ^l, T _y ^m), (T _x ^m, T _y ^h), { (T _x ^m, T _y ^m), (T _x ^h, T _y ^h);

C) judge whether R comprises target translation transfer function.Judgment task carries out in two steps: at first calculate in unit R to be selected each marginal point P ∈ E ₂Pairing indefinite zone; If the zone comprises at least one point of interest (satisfying the marginal point of relational expression (2)), then this zone is decided to be qualified zone; Secondly, verify in unit R, qualified number of regions accounts for the ratio of total edge point; If this ratio value is greater than predefined fth value, unit R then to be selected is designated as effectively.

D) to the effective subelement R of each layer _i, implement ordering according to the fth value, choose maximum and inferior big person and be pressed into storehouse, as the unit to be selected of further division.

E) iterative step a)-d), the length until indefinite zone converges on predefined τ value.The iterative processing of step a)-d) is equivalent in essence and makes up a binary search tree.

Below the extraction of characteristic angle point is described.

Character such as the feature angle point is an image representative point feature, and it has in image relative orientation obviously, the geometric position is accurate.Can select for use classical Forstner operator that main and auxiliary image is carried out the feature angle point extracts, this operator is by calculating Robert ' the s gradient and the pixel (c of each pixel, r) be the gray scale covariance matrix of a window (as 5 * 5) at center, in image, seek have as far as possible little and near the point of the error ellipse of circle as unique point, its step comprises:

(1) calculates Robert ' the s gradient of each pixel

$\left.\begin{array}{c}{g}_u=\frac{\∂g}{\∂u}=g_{i+1,j+1}-g_{i,j}\\ g_v=\frac{\∂g}{\∂v}=g_{i,j+1}-g_{i+1,j}\end{array}\right\}$

(2) covariance matrix of gray scale in calculating l * l (as 5 * 5 or bigger) window

$Q=N^{-1}={\left[\begin{array}{cc}\mathrm{\&Sigma;}{{\mathrm{<figure-callout\; id="2"\; label="g\; u"\; filenames="A20051013220000221.png"\; state="\{\{state\}\}">g</figure-callout>}}_u}^2& \mathrm{\&Sigma;}{g}_u{g}_v\\ \mathrm{\&Sigma;}{g}_v{g}_u& \mathrm{\&Sigma;}{g_v}^2\end{array}\right]}^{-1}$

Wherein

$\mathrm{\&Sigma;}{{\mathrm{<figure-callout\; id="2"\; label="g\; u"\; filenames="A20051013220000221.png"\; state="\{\{state\}\}">g</figure-callout>}}_u}^2=\underset{i=c-k}{\overset{c+k-1}{\mathrm{\&Sigma;}}}\underset{j=r-k}{\overset{r+k-1}{\mathrm{\&Sigma;}}}{(g_{i+1,j+1}-g_{i,j})}^2$

$\mathrm{\&Sigma;}{{\mathrm{<figure-callout\; id="2"\; label="g\; v"\; filenames="A20051013220000221.png"\; state="\{\{state\}\}">g</figure-callout>}}_v}^2=\underset{i=c-k}{\overset{c+k-1}{\mathrm{\&Sigma;}}}\underset{j=r-k}{\overset{r+k-1}{\mathrm{\&Sigma;}}}{(g_{i,j+1}-g_{i+1,j})}^2$

$\mathrm{\&Sigma;}{g}_u{g}_v=\mathrm{\&Sigma;}{g}_v{g}_u=\underset{i=c-k}{\overset{c+k-1}{\mathrm{\&Sigma;}}}\underset{j=r-k}{\overset{r+k-1}{\mathrm{\&Sigma;}}}(g_{i+1,j+1}-g_{i,j})(g_{i,j+1}-g_{i+1,j})$

k＝INT(l/2)

(3) calculate interest value q and w:

$q=\frac{4\mathrm{DetN}}{{\left(\mathrm{trN}\right)}^2}$

$w=\frac{1}{\mathrm{trQ}}=\frac{\mathrm{DetN}}{\mathrm{trN}}$

Wherein DetN represents the determinant of matrix N; TrN represents the mark of matrix N.

(4) determine to treat reconnaissance.If point of interest is greater than given threshold value, then this pixel point is for treating reconnaissance.Threshold value is an empirical value, can be with reference to as follows:

$\left.\begin{array}{c}{T}_q=0.5-0.75\\ T_w=\left\{\begin{array}{cc}f\stackrel{\&OverBar;}{w}& (f=0.5-1.5)\\ {\mathrm{cw}}_c& (c=5)\end{array}\right.\end{array}\right\}$

Wherein Be the weight average value; w _cIntermediate value for power.As q＞T _qWhile w＞T _wThe time, this pixel is for treating reconnaissance.

(5) with the weight w be foundation, select extreme point, promptly in certain suitably window, select the reconnaissance for the treatment of of w maximum.

Fig. 4 is for the flow chart of the same place intersection search that adopted in the method for the present invention, with reference to Fig. 4.

Image feature is meant the point of discontinuity of image greyscale curved surface.Point feature and line feature are the common characteristic informations of image.The point feature mainly indicates apparent point, as angle point, round dot etc.The line feature is meant image " edge ".So-called feature angle point promptly is image two " edge " point of crossing; On remote sensing image, the general corresponding road junction of feature angle point, building angle point, point of crossing, river or the like.

Selecting for use classical Forstner operator that main and auxiliary image is carried out the feature angle point extracts, this operator is by Robert ' the s gradient of each pixel of calculating with pixel (c, r) be the gray scale covariance matrix of a window (as 5 * 5) at center, in image, seek have as far as possible little and near the point of the error ellipse of circle as unique point, detailed algorithm is referring to pertinent literature.In order to raise the efficiency, only the feature angle point to be carried out in the pairing main and auxiliary overlapping region of the thick coupling of feature and extract; In the leaching process, adopted even Distribution Strategy, and given corresponding overlay area call number for simultaneously each feature angle point, the equalization that even distribution can make the feature angle point reach on the global sense distributes.Evenly Distribution Strategy comprises: at first, and to the even layout rules graticule mesh of main and auxiliary image overlap area; Secondly, use characteristic point extracts operator, extracts a feature angle point in each regular subnet.The overlay area call number numerically is equivalent to regular sublattice net numbering.

So-called same place promptly is the unique point of corresponding same atural object on the main and auxiliary image; The same characters of ground object point of main and auxiliary image is right to constituting same place, and in based on image coupling document, same place is to also can be described as the reference mark sometimes to (Control Points or Tie-Points).The sequential similarity measure method of provincial characteristics coupling is adopted in the same place search.Because slightly mate having carried out global characteristics before, the same place hunting zone dwindles greatly, and then bring descending significantly and mating the raising of correctness of computation complexity.In order to enrich same place information, adopt main and auxiliary image intersection search method, may further comprise the steps:

With main image feature angle point is reference, according to the transformation relation of the thick coupling of feature, seeks same place in the auxilliary image field of search; Then the gained same place to giving corresponding call number, i.e. same place image coordinate residing overlay area call number in main and auxiliary image.

With auxilliary image feature angle point is reference, a same step cross processing.

The gained same place is combined processing.In the merging process, two groups of same places are identical to call number, then only get one; Otherwise both all are selected in, and include final same place in in the tabulation.

The employing of intersection search method has been considered main and auxiliary image character information comprehensively with waiting power, and it is right further to have enriched same place, for the enforcement based on the adaptive region block algorithm of least square is got ready.

For the high resolution SAR data, shade is a big key character of image.Classical feature angle point extraction algorithm can not well be removed shadow character point (shadow character is put non-rigid ground object target, does not satisfy the requirement as same place), and this obviously can influence the right accuracy of overall same place to a certain extent.Therefore, remove high resolution SAR image same place centering shadow effect and will seem especially important.In addition, this method also adopts the thought of probability statistics, to same place to having carried out the elimination of rough difference processing, thereby further guaranteed finally to obtain the right accuracy of same place.

Fig. 5 is the flow chart of the LS-MMSE algorithm flow that adopted in the method for the present invention.

Based on least square region adaptivity block algorithm is the gordian technique of guaranteeing the sub-pixel level precision index.High resolution SAR image map sheet coverage is big, and SAR image low coverage contraction that local different incidence angles is introduced and sensing system error etc. all are to cause that the direct error source of local deformation takes place image.Therefore, auxilliary image integral body is carried out single geometric transformation and will become unrealistic.

Auxilliary image geometry rectification of distortion method adopts the affined transformation method, sees formula (III):

x ₂＝a ₀+a ₁x+a ₂y

(III)

y ₂=b ₀+b ₁x+b ₂y

Why select affined transformation, this has its characteristic decision.At first, this conversion simple, intuitive; Secondly, affined transformation can not only be used for handling airborne image, and is applicable to spaceborne image (polynomial expression of high-order is corrected and generally only is applicable to spaceborne data); At last, the comprehensive stack of piecemeal subimage local affine transformations on image overall situation yardstick, is equivalent to a high-order geometric transformation function of match formula.

Concrete steps based on least square region adaptivity block algorithm are as follows:

1) prime area piecemeal uses formula (III), obtains initial piecemeal subimage same place to error equation (IV):

$V_x^i=\frac{\&PartialD;x}{\&PartialD;a_0^i}{\mathrm{da}}_0^i+\frac{\&PartialD;x}{\&PartialD;a_1^i}{\mathrm{da}}_1^i+\frac{\&PartialD;x}{\&PartialD;a_2^i}{\mathrm{da}}_2^i-l_x^i$

$V_y^i=\frac{\&PartialD;y}{\&PartialD;b_0^i}{\mathrm{db}}_0^i+\frac{\&PartialD;y}{\&PartialD;b_1^i}+\frac{\&PartialD;y}{\&PartialD;b_2^i}{\mathrm{db}}_2^i-l_y^i---\left(\mathrm{IV}\right)$

$l_x^i=x_2^i-(a_0^i+a_1^{i}x+a_2^{i}y)$

$l_y^i=y_2^i-(b_0^i+b_1^{i}x+b_2^{i}y)$

Formula (IV) is expressed as matrix form, obtains formula (V):

V=AX-L

In the formula $V=\left[\begin{array}{c}{V}_x^i\\ V_y^i\end{array}\right],$ $A=\left[\begin{array}{cccccc}1& x& y& 0& 0& 0\\ 0& 0& 0& 1& x& y\end{array}\right],$ $X=\left[\begin{array}{c}{\mathrm{da}}_0^i\\ {\mathrm{da}}_1^i\\ {\mathrm{da}}_2^i\\ {\mathrm{db}}_0^i\\ {\mathrm{db}}_1^i\\ {\mathrm{db}}_2^i\end{array}\right],$ $L=\left[\begin{array}{c}{l}_x^i\\ l_y^i\end{array}\right].$

(V)

2) method error equation (V):

(A ^TA)X＝A ^TL

X＝(A ^TA) ^-1A ^TL

Iterative computation deformation parameter corrected value da ₀ ⁱ, da ₁ ⁱ, da ₂ ⁱ, db ₀ ⁱ, db ₁ ⁱ, db ₂ ⁱAnd corresponding deformation parameter a ₀ ⁱ, a ₁ ⁱ, a ₂ ⁱ, b ₀ ⁱ, b ₁ ⁱ, b ₂ ⁱThrough deriving current deformation parameter a ₀ ⁱ, a ₁ ⁱ, a ₂ ⁱ, b ₀ ⁱ, b ₁ ⁱ, b ₂ ⁱCan be by a preceding deformation parameter a ₀ ^I-1, a ₁ ^I-1, a ₂ ^I-1, b ₀ ^I-1, b ₁ ^I-1, b ₂ ^I-2Calculate acquisition via relational expression (VI):

$a_0^i=a_0^{i-1}+{\mathrm{da}}_0^i+a_0^{i-1}{\mathrm{da}}_1^i+b_0^{i-1}{\mathrm{da}}_2^i$

$a_1^i=a_1^{i-1}+a_1^{i-1}{\mathrm{da}}_1^i+b_1^{i-1}{\mathrm{da}}_2^i$

$a_2^i=a_2^{i-1}+a_2^{i-1}{\mathrm{da}}_1^i+b_2^{i-1}{\mathrm{da}}_2^i$

$b_0^i=b_0^{i-1}+{\mathrm{db}}_0^{i-1}+a_0^{i-1}{\mathrm{db}}_1^i+b_0^{i-1}{\mathrm{db}}_2^i---\left(\mathrm{VI}\right)$

$b_1^i=b_1^{i-1}+a_1^{i-1}{\mathrm{db}}_1^i+b_1^{i-1}{\mathrm{db}}_2^i$

$b_2^i=b_2^{i-1}+a_2^{i-1}{\mathrm{db}}_1^i+b_2^{i-1}{\mathrm{db}}_2^i$

3) subimage block division, rollback increases the reference mark and handles.

Division: when step 2) deformation parameter is separated when being tending towards restraining in, calculates the l of this subimage block correspondence _x ⁱ, l _y ⁱ, if when calculating institute's value greater than assign thresholds (sub-pixel level precision, this threshold value are 1), be the dad image piece with this image block then, one is divided into four, division generates four littler number of sub images pieces of overlay area.

Rollback: behind the dad image block splitting, newly obtain same place that subimage block comprises and during less than necessity observation number, give the deformation parameter value of dad image piece correspondence then directly for this new subimage block to (reference mark) number.Level and smooth for the image subimage block edge transition that guarantees to resample, the Algorithm Error precision reaches MMSE (Minimum Mean Square Error), often needs the excess observation of some." rollback " process is harmful to final matching precision.Its essence is hinting that predetermined matching precision is high more, and then the number of control points of demand is also many more.Therefore, to count be the most direct, the effective means of avoiding " rollback " effect in Chong Zu control.

Increase the reference mark: in order to avoid the rollback effect as far as possible, when same place in the subimage block is observed number to (reference mark) number less than necessity, again extract the feature angle point and encrypt same place, and then can continue the whole iterative processing in subimage block division, rollback and increase reference mark (reference mark).It should be noted that to comprise quantity of information limited in view of subimage block, increase the reference mark and handle and often only need iteration 1～2 time; This means that also for the poor relatively zone of image information, the rollback effect is inevitable.

4) revise the right image coordinate value of same place.When the corresponding subimage block division of auxilliary image, rollback and increase reference mark process are ended, interative computation step 1) again, 2); But, this moment need be to x ₂ ⁱ, y ₂ ⁱDo following correcting process, until all same places to disposing:

$x_2^i=x_2^{i-1}+V_x^{i-1}$

$y_2^i=y_2^{i-1}+V_y^{i-1}$

5) adopt indirect method, use corresponding deformation parameter to carry out geometric transformation to the subimage block that comprises in the auxilliary image successively, finally obtain the image (resample and select bilinear interpolation) behind the registration.

Method synthesis of the present invention has been considered the robustness of characteristic matching and the accuracy of zone coupling, the overall employing: based on the thick matching algorithm of Hausdorff distance measure global characteristics of binary tree search by slightly implementing autoregistration image to smart matching strategy, solved the counting yield bottleneck problem of Hausdorff distance measure well, algorithm is sane, reached the matching precision of sub-pixel-level, thereby dwindled the indefinite zone of same place search, improved treatment effeciency.The same place coupling has adopted sequential similarity to detect; The introducing of intersection search method waits power to consider the characteristic information of main and auxiliary image, and it is right to have enriched same place; The same place coupling has obtained the matching precision of Pixel-level between image.At the same place that obtains the Pixel-level precision under the prerequisite, proposition based on least square region adaptivity block algorithm, be to guarantee to reach sub-pixel-level on the image global sense (to adopt the 3 meters carried SAR images in 9 air strips, area, Xinghua, Jiangsu in July, 2003 as experimental data, registration accuracy is (0.810076,0.833336) pixel) key link of matching precision; The enforcement of adaptive region partition strategy has further increased the dirigibility and the robustness of algorithm.

High resolution SAR image registration processing method of the present invention is specifically used

Application example one:

Data constitute: area, Xinghua, Jiangsu in July, 2,003 9 air strips carried SAR images (CAS Electronics Research Institute's carried SAR data, resolution 3m).

Processing platform: China Remote Satellite Earth Station, CAS's microwave group is ground the CAESAR-RegSAR system software voluntarily.

Effect: realized 9 the automatic/semi-automatic registrations of air strips Airborne High-resolution SAR image in area, Xinghua, Jiangsu, generated big map sheet and inlay image.

Effect: reached overall sub-pixel-level registration accuracy, registration accuracy is (0.810076,0.833336) pixel.

Application example two:

Data constitute: area, Sanya, Hainan in January, 2005 Ku, X, L, 5 band carried SARs of P four wave bands image (China Electronic Science and Technology Corporation's the 38th research institute's carried SAR data, resolution 1m).

Processing platform: the CAESAR-RegSAR system software that China Remote Satellite Earth Station, CAS's microwave group is researched and developed voluntarily.

Effect: realized Hainan three subregion Ku, X, L, the automatic/semi-automatic registration of P four wave bands 5 band carried SARs image respectively, generated big map sheet and inlay image.

Effect: reached overall sub-pixel-level registration accuracy, registration accuracy is (0.730816,0.582958) pixel.

Claims (7)

1, a kind of high resolution SAR image registration processing method is characterized in that, comprises step:

(a) create pyramid image, extract the image texture feature, obtain marginal information as coupling input data;

(b) carry out the thick coupling of the image overall situation based on the image texture feature of being extracted, obtain the main and auxiliary image behind the thick registration;

(c) the feature angle point of main and auxiliary image behind the thick registration of extraction, and carry out the same place intersection search, it is right to obtain same place;

(d) handle carrying out the adaptive region piecemeal according to same place.

2, high resolution SAR image registration processing method as claimed in claim 1 is characterized in that, wherein step (a) comprising: create pyramid image; Adopt the Canny algorithm, main and auxiliary image is carried out edge extracting; Edge image is removed noise processed; According to edge length described edge image is further handled.

3, high resolution SAR image registration processing method as claimed in claim 1 is characterized in that, wherein adopts quick Hausdorff distance algorithm slightly to mate in the step (b).

4, high resolution SAR image registration processing method as claimed in claim 1 is characterized in that, wherein adopts even Distribution Strategy to extract the feature angle point in the step (c), comprising: to the even layout rules graticule mesh of main and auxiliary image overlap area; Use characteristic point extracts operator, extracts a feature angle point in each regular subnet.

5, high resolution SAR image registration processing method as claimed in claim 1, it is characterized in that, wherein adopt main and auxiliary image intersection search method in the step (c), be reference with main and auxiliary image feature angle point respectively, according to the transformation relation of the thick coupling of feature, in the image field of search, seek same place; Then the gained same place to giving corresponding call number; With the gained same place to merging processing.

6, high resolution SAR image registration processing method as claimed in claim 1 is characterized in that, wherein comprises step in step (d):

(d1) prime area piecemeal obtains initial piecemeal subimage same place to error equation (IV) according to formula (III),

x ₂＝a ₀+a ₁x+a ₂y (III)

y ₂＝b ₀+b ₁x+b ₂y

$V_x^i=\frac{\&PartialD;x}{{\&PartialD;a}_0^i}{\mathrm{da}}_0^i+\frac{\&PartialD;x}{{\&PartialD;a}_1^i}{\mathrm{da}}_1^i+\frac{\&PartialD;x}{{\&PartialD;a}_2^i}{\mathrm{da}}_2^i-l_x^i$

$V_y^i=\frac{\&PartialD;y}{{\&PartialD;b}_0^i}{\mathrm{db}}_0^i+\frac{\&PartialD;y}{{\&PartialD;b}_1^i}{\mathrm{db}}_1^i+\frac{\&PartialD;y}{{\&PartialD;b}_2^i}{\mathrm{db}}_2^i-l_y^i---\left(\mathrm{IV}\right)$

$l_x^i=x_2^i-(a_0^i+a_1^{i}x+a_2^{i}y)$

$l_y^i=y_2^i-(b_0^i+b_1^{i}x+b_2^{i}y)$

(d2) iterative computation deformation parameter corrected value da ₀ ⁱ, da ₁ ⁱ, da ₂ ⁱ, db ₀ ⁱ, db ₁ ⁱ, db ₂ ⁱAnd corresponding deformation parameter a ₀ ⁱ, a ₁ ⁱ, a ₂ ⁱ, b ₁ ⁱ, b ₁ ⁱ, b ₂ ⁱ

(d3) subimage block division, rollback increases the reference mark and handles;

(d4) revise the right image coordinate value of same place;

(d5) use corresponding deformation parameter to carry out geometric transformation to the subimage block that comprises in the auxilliary image successively, obtain the image behind the registration.

7, a kind of high resolution SAR image registration disposal system is characterized in that, comprising:

The image data extraction unit is used to extract the image texture characteristic, obtains marginal information and this information is sent to the image pretreatment unit as thick coupling input data;

The image pretreatment unit carries out the thick coupling of the overall situation to the image that extracts the image texture feature, obtains the main and auxiliary image information behind the registration, and relevant information is offered the Image registration processing unit processes;

The Image registration processing unit is used to extract the feature angle point of the main and auxiliary image that obtains after the thick coupling, and carries out the same place intersection search, and it is right to obtain same place; Handle carrying out the adaptive region piecemeal according to same place.

Images