CN106093939A - A kind of InSAR image phase unwrapping method based on phase contrast statistical model - Google Patents

Abstract

The present invention relates to Radar Technology field, particularly relate to the phase unwrapping technical field in synthetic aperture radar interferometry (InSAR) image procossing.The present invention constructs a kind of InSAR image phase unwrapping method based on phase contrast statistical model, first with phase contrast derivative function and Terrain Elevation variable quantity relational expression, deriving InSAR phase contrast statistical function model, system decoherence factor and the hypsography change impact on phase unwrapping are taken in the change of phase place into account to utilize phase contrast statistical function to estimate；In phase unwrapping, utilize phase contrast statistical model as the Nonlinear Constraints of least square phase unwrapping, and the extreme-value problem that solution twines in Models computed before and after's phase contrast replaces with solution and twines the integer multiple difference problem of before and after's phase contrast, carry out phase unwrapping by integer multiple difference minima before and after solving phase unwrapping.

Description

A kind of InSAR image phase unwrapping method based on phase contrast statistical model

Technical field

The present invention relates to Radar Technology field, particularly relate in synthetic aperture radar interferometry (InSAR) image procossing 's

Phase unwrapping method.

Background technology

Phase unwrapping is one of committed step of InSAR data process, and the precision of its data processed result directly influences The precision of final data product, but phase unwrapping is referred to as np problem, by two class phase place decoherence factors, (phase place is made an uproar for it Sound) impact: a class is that processing system itself causes, and such as thermal noise, Temporal decoherence, baseline decoherence are (due to antenna Sight line is the most parallel), and interpolation algorithm or Image filter arithmetic and the noise that introduces in process of image registration；Another kind of it is The fluctuations of landform causes phase place lack sampling to cause geometric phase distortion etc., causes coherent speckle noise.

The function model of the interferometric phase Φ obtained after being processed by two width radar images is write as linear term φ_linearAnd non-thread Property item φ_non-linearThe form of sum, such as formula (1):

$\begin{array}{c}\mathrm{\Φ}={\mathrm{\φ}}_{linear}+{\mathrm{\φ}}_{non-linear}=\frac{4\mathrm{\π}}{\mathrm{\λ}}{D}_{linear}+\frac{4\mathrm{\π}}{\mathrm{\λ}}{D}_{non-linear}\\ =\frac{4\mathrm{\π}}{\mathrm{\λ}}(\frac{B_{\⊥}}{R\sin \mathrm{\θ}}\mathrm{\Δ}h+v\·T)+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·D_{non-linear}\end{array}---\left(1\right)$

(1) in formula, D_linear、D_non-linearRepresent that hypsography changes the linear processes part caused, v table respectively Linear Ground Deformation speed, T represents the radar data time reference line interfered obtaining, and λ is radar wavelength, and R is that oblique distance is long, Δ H is that twice observed altitude of tested point is poor, B_⊥It is being perpendicular to the projection of direction of visual lines, i.e. effective base length for baseline.Formula (1) First two on the right of middle equal sign is linear segment, linear with vertical parallax, time reference line；Section 3 is non-linear Point, non-linear partial causes tradition branch tangent line solution to twine loop and cannot close and solution twines the phenomenons such as error propagation, directly with linearly Method cannot obtain high-precision disentanglement fruit.

According to whether take phase nonlinear into account, phase unwrapping method is divided into linear processes phase unwrapping method.At present, The most representational in linear phase unwrapping method is least square (Least Squares, LS) unwrapping method and weighting minimum Two take advantage of (Weighted Least Squares, WLS) phase unwrapping method.Refer to document: Lu Y G, Wang X ZH, Zhang X P.Weighted least-squares phase unwrapping algorithm based on derivative variance correlation map.2007.Linear method is based on a supposed premise: be wound around the discrete partial derivative of phase place (i.e. adjacent to the phase contrast of pixel) presents the linear trends of change of increasing or decreasing, and the absolute value of phase contrast is both less than π, to discrete Partial derivative just can try to achieve the phase place after solution twines along any direction iterated integral.Linear method can be reduced by phase weighting The phase masses that coherence factor causes reduces the impact twining solution, to a certain extent, it is to avoid error is in overall situation phase diagram Transmission, improves the speed of model convergence by the method for iteration.But, linear phase solution twines algorithm and does not differentiates between identification phase noise Source, be easily caused global error transmission, it is impossible to avoid the interference fringe caused due to orographic factor to lose or the folded feelings covered Condition, fractional phase can not get abundant solution to be twined, and the precision of solving result is the highest；Weighted phases solution twines and avoids merely phase residual error point Solution is twined the impact in region, it is impossible to avoiding the interference fringe caused due to orographic factor to lose or folded situation about covering, reason is shadow Ring the deviation essence that the decoherence factor such as phase place geometric distortion of phase masses, phase noise etc. make phase unwrapping result occur On be that solution twines model and do not accounts for the phase place of nonlinear change and the deviation that produces.

Currently for nonlinear optimization method, actually construct nonlinear model the process solved, the most non- Linear optimization estimation problem.Phase unwrapping based on optimization method is when being iterated solving, it is desirable to pre-estimates solution and twines phase place Initial value, then by the method acquisition algorithm of linear approximation or the optimal solution of model.Refer to technical literature 1:Wang Q S, Huang H F,An X Y,et al.Improving phase Unwrapping Techniques by the Use of Nonlinear Phase Model 2011,2:Eineder M, Adam N.A maximum likelihood estimator To simultaneously unwrap 2005, disclosed in technology contents.But, nonlinear solution twines model and is solving big picture During element interference image, use the alternative manner of linear approximation, it is impossible to improve efficiency and the precision of disentanglement fruit of computing, even have Shi Wufa solves, and in filtering or solution twine, the most directly takes the topography variation impact on phase contrast into account, and the result obtained is also Inaccuracy.

Summary of the invention

The present invention constructs InSAR image phase unwrapping method based on phase contrast statistical model, first with phase contrast Derivative function and Terrain Elevation variable quantity relational expression, derive InSAR phase contrast statistical function model, utilizes phase contrast to add up letter Number estimates that system decoherence factor and the hypsography change impact on phase unwrapping are taken in the change of phase place into account；At phase unwrapping In, utilize phase contrast statistical model as the Nonlinear Constraints of least square phase unwrapping, and in Models computed will solve Twine the extreme-value problem of before and after's phase contrast to replace with solution and twine the integer multiple difference problem of before and after's phase contrast, by solving phase unwrapping Front and back integer multiple difference minima carries out phase unwrapping.

Following concrete technical scheme is employed for realizing foregoing invention content:

Step 1, suppose that pending interference image row, column pixel size is M × N, be divided into the local of k × k size Estimate window, calculate phase place in difference DELTA φ of ranks both direction, and to phase contrast in orientation to distance to derivation, obtain Phase difference derivative (phase frequency under discrete phase)

$\frac{{\∂}^2\mathrm{\φ}}{\∂y\∂x}=\sqrt{\underset{i=1}{\overset{M-1}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{({\mathrm{\Δx}}_{i,j}-\underset{k\×k}{E}(\mathrm{\Δ}x\left)\right)}^2}+\sqrt{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N-1}{\mathrm{\Σ}}}{({\mathrm{\Δy}}_{i,j}-\underset{k\×k}{E}(\mathrm{\Δ}y\left)\right)}^2}/\left(k^2\right),(i=1,...M-1;j=1,...N-1);$

Wherein, Δ x_i,j=x_i+1,j+x_i-1,j-2x_i,j, Δ y_i,j=y_i,j+1+y_i,j-1-2y_i,jRepresent that interferogram is along x, y respectively The phase difference value in direction；It is the phase contrast average in x, y direction in k × k window respectively；It is with window The phase contrast derivative standard deviation value at center.

Step 2, according to upper step phase contrast derivative value, calculate phase contrast statistical model in partial estimation window, and estimation office The distribution of Terrain Elevation variable quantity in portion's window；

Step 3, relational expression according to phase contrast frequency and terrain slope angle estimate Terrain Elevation variable quantity；

Phase contrast statistical model in step 4, acquisition partial estimation window；

Step 5, estimating window intra-orally shape variable quantity distribution function is utilized phase contrast in window to be smoothed, at estimating window In Kou, distance is multiplied by phase contrast statistical model to orientation to phase difference value, in available window, takes the phase of topography variation amount into account Potential difference estimates model.

Step 6, solution is twined longitudinal separation to and orientation to integer multiple difference DELTA k of phase contrast_yWith Δ k_xCarry out estimating also Integration, it is possible to obtain the phase contrast integer after solution twines, is converted to solve integer multiple difference by the problem solving phase contrast Problem, the functional equation of structure integer multiple parameter is:With Wherein Δ φ_yWith Δ φ_xRepresent that phase point is wound around phase contrast,WithRepresent that solution twines rear phase contrast；

Step 7, least square phase unwrapping restricted model is replaced with integer multiple difference constrained minimization model, such as public affairs Formula:

$\min \left(J\right)=\{\underset{i=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}{\left(P\right({\mathrm{\Δ\φ}}_x\left)\mathrm{\Δ}k\right)}_{i,j}^2+\underset{i=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}{\left(P\right({\mathrm{\Δ\φ}}_y\left)\mathrm{\Δ}k\right)}_{i,j}^2\}$

$\begin{array}{ccc}s.t& \underset{K}{\&Sigma;}\&lsqb;\left(\right({\mathrm{Pk}}_{i,j+1}^y-{\mathrm{Pk}}_{i,j}^y)-({\mathrm{Pk}}_{i+1,j}^x-{\mathrm{Pk}}_{i,j}^x))({\mathrm{\&Delta;k}}^t-{\mathrm{\&Delta;k}}^{t-1})\&rsqb;<\mathrm{\&xi;}& t>1,(i,j)\&Element;S_{k\&times;k}\end{array}$

Wherein, P (Δ φ_x)、P(Δφ_y) represent phase point row vector, column vector phase contrast statistical value, Δ k^tWith Δ k^t-1 Representing that the integer multiple of twice iteration is poor, ξ is minimum；Constraints s.t represents in local window, the t time and the t-1 time Solution during iterative twines integer multiple difference and meets the requirement that difference is minimum, and K represents the size estimating window；

Step 8, utilize discrete least-squares iteration method obtain solution twine phase place.

The beneficial effect that technical solution of the present invention is brought

(1) InSAR phase unwrapping method based on phase contrast statistical model, utilizes interference phase difference statistical model to take ground into account Shape factor, can play the purpose of pre-estimation in solution twines to topography variation.For alpine terrain and breaks landform with And when landform phase place has sedimentation, for the interference fringe of formation comprises the interference data of aliasing phase place, add up based on phase contrast The noise that feature detection landform causes, the solution of acquisition twines image smoothing, and weight Twined Structure is more nearly original interferogram, and solution twines essence Du Genggao.

(2) present invention uses phase contrast statistical value to replace phase difference value as the additional constraint bar of least square phase unwrapping Part, phase contrast affected by noise can be screened by additional constraint condition, eliminates noise while ensureing smoothing pseudorange Solution is twined the impact of path of integration.

(3) in solution twines model solution, for processing complicated Hessian matrix and the problem of big matrix convergence, and Take into account and how to change constraints, twine effect in the hope of reaching optimum solution under conditions of precision ensureing, solve relatively difficult, The problem that effect is undesirable, the present invention is by public for the least square phase unwrapping equivalent model that model conversion is additional constraint condition Formula, is converted to the problem solving phase difference estimation solve the problem that integer multiple is estimated, utilizes phase contrast statistical eigenfunction Estimating the result of each iteration, iteration all embodies discrete optimized step each time.

Accompanying drawing explanation

Fig. 1 is that different incidence angles affects figure to the identical gradient；

Fig. 2 is that InSAR interference phase difference statistical nature model obtains flow process.

Marginal data, represents different side looking radar incident angles, θ 1., 2., the most respectively in Fig. 1₁、θ₂、θ₃Represent respectively Survey terrain slope inclination angle, district；Shadow is the shadow region without radar appearance.

Detailed description of the invention

Below in conjunction with the accompanying drawings technical solution of the present invention is further illustrated.

In the present invention, technology describes for convenience, has first carried out defined below:

Definition 1: set interferometer radar and observe oblique distance difference for twice as Δ R, twice observation phase contrast is Δ φ formula, both public affairs Formula is as follows,

$\mathrm{\&Delta;}R=B_{//}+\frac{B_{\&perp;}}{R_0}\&lsqb;\frac{R_1-R_0}{tan({\mathrm{\&theta;}}_0-{\mathrm{\&tau;}}_y)}+\frac{(x-x_0){\mathrm{tan\&tau;}}_x+h}{\sin ({\mathrm{\&theta;}}_0-{\mathrm{\&tau;}}_y)}{\mathrm{cos\&tau;}}_y\&rsqb;---(1-1)$

$\mathrm{\&Delta;}\mathrm{\&phi;}=\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\{\frac{B_{\&perp;}}{R_0}\frac{1}{tan({\mathrm{\&theta;}}_0-{\mathrm{\&tau;}}_y)}(R_1-R_0)+\frac{B_{\&perp;}}{R_0}\frac{{\mathrm{tan\&tau;}}_x{\mathrm{cos\&tau;}}_y}{sin({\mathrm{\&theta;}}_0-{\mathrm{\&tau;}}_y)}(x-x_0)+\frac{B_{\&perp;}}{R_0}\frac{{\mathrm{cos\&tau;}}_y}{sin({\mathrm{\&theta;}}_0-{\mathrm{\&tau;}}_y)}h\}---(1-2)$

Wherein, λ is radar wavelength, R₀、R₁Being respectively twice oblique distance long, h is tested point height；B_⊥It is being perpendicular to for baseline The projection of direction of visual lines, i.e. effective base length；B_//For horizontal base line；(x₀,y₀,z₀) and (x, y z) are respectively twice observation Coordinate；τ_x、τ_yFor ground observation point along orientation to distance to the angle of gradient, θ₀It it is the meansigma methods of twice observation radar downwards angle of visibility.

The most in fig. 1: in the region that gradient fluctuations is bigger, angle of incidence is more than θ₂Time radar image will produce the moon Shadow, along with angle of incidence reduces, shadow region is also being shunk；Little to θ at angle of incidence₃Time will appear from push up the inverted phenomenon in the end cause radar The geometric distortion of image.The terrain slope different when angle of incidence is identical is also identical on the impact of phase contrast.

Definition 2: Terrain Elevation variable quantity formula

Hypsography change for observation station can utilize its point high variable quantity to estimate, and Terrain Elevation becomes Change amount is it is believed that its logarithmic function obeying three parameter impacts is distributed:

$\log L_{\mathrm{\&alpha;},v,\mathrm{\&mu;}}\left(\right|d\left|\right)=\frac{L_{\mathrm{\&alpha;},v,\mathrm{\&mu;}}\&lsqb;log\left(\right|d\left|\right)\&rsqb;}{\left|d\right|}---(1-3)$

Wherein, | d | represent Terrain Elevation variable quantity d projection from radar line of sight to direction (side looking radar to sight line to Phase place change most sensitive): d=| d | cos θ (θ represent observation station tangent line and radar line of sight to angle), | d |=[tan τ_ytanτ_x]^T=[d_yd_x]^T.α, ν, μ are that three functions control parameter, can be right by changing three parameters for different landform Terrain Elevation change carries out preferable matching.

The explanation of parameter and value be can be found in following table:

Concrete technical scheme steps is as follows:

Step 1, the preparation of basic data.

Obtain two width SAR image registrations, filtered interference image, calculate coherence factor γ and discrete interferometric phase φ, Including the system parameter message (if wavelength, baseline, downwards angle of visibility, oblique distance and the landform DEM information that can obtain) added.In reality In the image processing process of border, common-used formulaCarry out interfering the estimation of coherent value,Represent SAR image S of two width areals after filtering₁With S₂Conjugate multiplication.In formula, l represents that regard counts more.And interferogram is read Being taken as the phasing matrix of relocatable, and labelling ranks size is M, N, unit is pixel.

Phase contrast derivative value in step 2, calculation interferogram row, column direction partial estimation window.

Interference image is divided into the estimation window of k × k, calculates phase place difference DELTA φ in ranks both direction, and right Phase contrast in orientation to distance to derivation, obtain phase contrast derivative (the phase contrast frequency under discrete phase).At discrete picture In, phase derivative available phases difference represents:

$\frac{{\&part;}^2\mathrm{\&phi;}}{\&part;y\&part;x}=\sqrt{\underset{i=1}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\mathrm{\&Delta;x}}_{i,j}-\underset{k\&times;k}{E}(\mathrm{\&Delta;}x\left)\right)}^2}+\sqrt{\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N-1}{\mathrm{\&Sigma;}}}{({\mathrm{\&Delta;y}}_{i,j}-\underset{k\&times;k}{E}(\mathrm{\&Delta;}y\left)\right)}^2}/\left(k^2\right),(i=1,...M-1;j=1,...N-1)$

Wherein, Δ x_i,j=x_i+1,j+x_i-1,j-2x_i,j, Δ y_i,j=y_i,j+1+y_i,j-1-2y_i,jRepresent that interferogram is along ranks respectively The phase difference value in direction；It is the phase contrast average in ranks direction in k × k window respectively；It is with window The phase contrast derivative standard deviation value at mouth center.

Step 3, the phase contrast derivative value drawn according to above-mentioned steps 2, calculate phase contrast statistics mould in partial estimation window Type, and estimate the distribution of Terrain Elevation variable quantity in local window.

The distribution of mapping district Terrain Elevation variable quantity is obtained according to landform high variable quantity formula (1-3).To by three ginsengs The Terrain Elevation value variable quantity d that digital-to-analogue is intended takes the logarithm and removes their Jacobian, if the scope district of phase estimation window Between size be k × k, d_x、d_yFor Terrain Elevation variable quantity d distance to orientation to the projection of both direction, then partial estimation In window, the two-dimensional probability density function formula of Terrain Elevation variable quantity is:

${\mathrm{PDF}}_{d_x,d_y}(d_x,d_y)=\frac{\log L_{\mathrm{\&alpha;},n,\mathrm{\&mu;}}\left(\sqrt{{d_x}^2+{d_y}^2}\right)}{2\mathrm{\&pi;}J}---\left(2\right)$

In above formula, J is Jacobian, wherein

By estimating window intra-orally shape high variable quantity density integral can be obtained its distribution function, ifSo discrete distribution function F of interferometric phase image mesorelief high variable quantity_dLetter for variable Δ t Yu d Number, formula is as follows:

$F_d\left(d\right)=\underset{k\&times;k}{\&Sigma;}\frac{\log L_{\mathrm{\&alpha;},v,\mathrm{\&mu;}}\left(\mathrm{\&Delta;}t\right)}{\mathrm{\&pi;}\sqrt{{\mathrm{\&Delta;t}}^2-d^2}}---\left(3\right)$

Step 4, relational expression according to phase contrast derivative and terrain slope angle estimate Terrain Elevation variable quantity.

A variable d describing Terrain Elevation is there is, due to Terrain Elevation in above formula topography variation amount distribution function Change is directly proportional relation with the change at terrain slope angle, then can be by the relation of phase contrast derivative Yu terrain slope angle Formula represents this variable.I.e. terrain slope change and the relation of phase derivative, such as following formula (4), (5):

$\begin{array}{c}{\mathrm{\&Delta;f}}_y=\frac{\&part;\mathrm{\&phi;}}{\&part;y}=k_B\frac{1}{tan(\mathrm{\&theta;}-{\mathrm{\&tau;}}_y)}\\ =k_B\frac{1+{\mathrm{tan\&theta;tan\&tau;}}_y}{\tan \mathrm{\&theta;}-{\mathrm{tan\&tau;}}_y}\end{array}---\left(4\right)$

$\begin{array}{c}{\mathrm{\&Delta;f}}_x=\frac{\&part;\mathrm{\&phi;}}{\&part;x}=k_B\frac{{\mathrm{tan\&tau;}}_x{\mathrm{cos\&tau;}}_y}{sin(\mathrm{\&theta;}-{\mathrm{\&tau;}}_y)}\\ =k_B\frac{{\mathrm{tan\&tau;}}_x}{\sin \mathrm{\&theta;}-{\mathrm{cos\&theta;tan\&tau;}}_y}\end{array}---\left(5\right)$

Wherein,X, y represent respectively the distance of interferometric phase to orientation to, Δ f_x、Δf_yRepresent respectively The phase contrast derivative of both direction；C is radar wave propagation speed；B_⊥It is being perpendicular to the projection of direction of visual lines, the most effectively for baseline The length of base；λ is wavelength；R is oblique distance；θ is angle of incidence；τ is the angle of gradient of surface sample point instantaneously.(4), (5) formula are carried out Inverse transformation obtains formula (6):

$\left\{\begin{array}{c}{d}_x({\mathrm{\&Delta;f}}_x,W)={\mathrm{\&Delta;f}}_{x}W\\ d_y\left({\mathrm{\&Delta;f}}_y\right)=\frac{{\mathrm{\&Delta;f}}_y\tan \mathrm{\&theta;}-k_B}{{\mathrm{\&Delta;f}}_y+k_B\tan \mathrm{\&theta;}}={\mathrm{Wk}}_B/\cos \mathrm{\&theta;}+\tan \mathrm{\&theta;}\end{array}\right.---\left(6\right)$

Wherein, W=cos θ (d_y-tanθ)/k_B, above formula explanation phase contrast derivative is with phase frequency, terrain slope angle and to enter Firing angle is the function of parameter.Except the component d of Terrain Elevation variable quantity in formula_x、d_yOutside two unknown quantitys, remaining parameter is System known parameters, or the variate-value asked for by step above.

Phase contrast statistical model in step 5, acquisition partial estimation window.

By step 3 to survey the two-dimensional probability density function formula (2) of district Terrain Elevation variable quantity and terrain slope angle with The relation function formula (6) of phase contrast derivative, derives the phase contrast statistical model such as following formula of both direction:

$\begin{array}{c}{P}_y\left({\mathrm{\&Delta;f}}_y\right)=\underset{k\&times;k}{l}\&CenterDot;d_y\left({\mathrm{\&Delta;f}}_y\right)F_d\left(d_y\right)\&CenterDot;{\mathrm{Tc}}_{\mathrm{\&theta;}}\\ =\underset{k\&times;k}{l}\&CenterDot;\frac{{\mathrm{\&Delta;f}}_y\tan \mathrm{\&theta;}/k_B-1}{({\mathrm{\&Delta;f}}_y/k_B+\tan \mathrm{\&theta;})}{F}_d\left(d_y\right)\&CenterDot;\left|\sin (\mathrm{\&theta;}-\arctan d_y)\right|\end{array}---\left(7\right)$

$\begin{array}{c}{P}_x\left({\mathrm{\&Delta;f}}_x\right)=\underset{k\&times;k}{l}\&CenterDot;d_x({\mathrm{\&Delta;f}}_x,W)F_d\left(d_x\right)\&CenterDot;{\mathrm{Tc}}_{\mathrm{\&theta;}}\\ =\underset{k\&times;k}{l}\&CenterDot;\frac{{\mathrm{\&Delta;f}}_x\cos \mathrm{\&theta;}(d_y-\tan \mathrm{\&theta;})}{k_B}\&CenterDot;\left|\sin (\mathrm{\&theta;}-\arctan d_y)\right|\end{array}---\left(8\right)$

Wherein, l is proportionality constant, and by estimating that window size is affected in estimating phase model, actual estimated window is more Little impact is the least, can will estimate that window is set to 3 × 3 sizes.Tc_θFor estimating to sample in window ratio, this be one by incidence angle θ With terrain slope angle arctand_yThe amount simultaneously controlled, both differences the biggest explanation interference effect is the poorest needs raising sample rate, Use the trigonometric function sin (θ-arctand of both change of reflection_y) represent.

Step 6, utilize partial estimation window in topography variation amount distribution function phase contrast in window is smoothed.

According to phase contrast statistical function formula (7), (8), in estimating window, distance is multiplied by phase difference value to orientation Phase contrast statistical model takes the phase difference estimation model of topography variation amount into account in can get window, such as (9), (10):

$\underset{k\&times;k}{E}\left({\mathrm{\&Delta;\&phi;}}_y\right)=P{\left({\mathrm{\&Delta;\&phi;}}_y\right)}_{k\&times;k}\&CenterDot;(\frac{1}{2\mathrm{\&pi;}}\&CenterDot;\frac{{\mathrm{\&Delta;\&phi;}}_y}{{\mathrm{\&Delta;f}}_{y_{i,j}}})\frac{1}{2{\mathrm{\&pi;\&Delta;f}}_{y_{i,j}}},(i=0,1,...k-2;j=0,1,...k-1)---\left(9\right)$

$\underset{k\&times;k}{E}\left({\mathrm{\&Delta;\&phi;}}_x\right)=P{\left({\mathrm{\&Delta;\&phi;}}_x\right)}_{k\&times;k}\&CenterDot;(\frac{1}{2\mathrm{\&pi;}}\&CenterDot;\frac{{\mathrm{\&Delta;\&phi;}}_x}{{\mathrm{\&Delta;f}}_{x_{i,j}}})\frac{1}{2{\mathrm{\&pi;\&Delta;f}}_{x_{i,j}}},(i=0,1,...k-2;j=0,1,...k-1)---\left(10\right)$

In formulaRepresent and estimate the estimated value of phase contrast in window k × k, and the relation between them meets following formula Condition:

$\begin{array}{ccc}P\left(\right|\underset{k\&times;k}{E}\left({\mathrm{\&Delta;\&phi;}}_{x_{i,j},y_{i,j}}\right)|>\mathrm{\&pi;})& or& P(-\mathrm{\&pi;}<\underset{k\&times;k}{E}({\mathrm{\&Delta;\&phi;}}_{x_{i,j},y_{i,j}})<\mathrm{\&pi;})\end{array}---\left(11\right)$

Transformation relation from phase contrast statistical nature Yu phase difference value, it is possible to use phase contrast statistical nature estimates phase The variation tendency of potential difference, identifies the phase residual error point caused due to landform.

Step 7, structure InSAR image least square phase unwrapping model based on phase contrast statistical model.

The quadratic sum minimum of discrete phase potential difference before and after phase unwrapping is utilized to build least square phase unwrapping model, the most such as Formula (12):

$\min \left(J\right)=\underset{i=1}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\mathrm{\&psi;}}_{i+1,j}-{\mathrm{\&psi;}}_{i,j}-{\mathrm{\&Delta;\&phi;}}_{i,j}^x)}^2+\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N-1}{\mathrm{\&Sigma;}}}{({\mathrm{\&psi;}}_{i,j+1}-{\mathrm{\&psi;}}_{i,j}-{\mathrm{\&Delta;\&phi;}}_{i,j}^y)}^2$

$\begin{array}{c}\begin{array}{cc}s.t& \underset{K}{\mathrm{\&Sigma;}}\&lsqb;M_{\mathrm{\&Delta;}\mathrm{\&psi;}|x_{i,j},y_{i,j}}\left(\mathrm{\&Delta;}{\stackrel{\&OverBar;}{\mathrm{\&psi;}}}^t\right|x_{i,j},y_{i,j})-M_{\mathrm{\&Delta;}\mathrm{\&psi;}|x_{i,j},y_{i,j}}\left(\mathrm{\&Delta;}{\stackrel{\&OverBar;}{\mathrm{\&psi;}}}^{t-1}\right|x_{i,j},y_{i,j})\&rsqb;<\mathrm{\&xi;}\end{array}\\ (i=0,1,...M-2;j=0,1,...N-1)\end{array}---\left(12\right)$

Wherein, M and N is respectively columns and the line number of interferogram phasing matrix, Δ φ_i,jIt is phase point (i, j) winding at place Phase contrast, ψ_i,jIt is that (i, j) solution at place twines rear phase place to phase point；ξ is minimum；Represent estimating window The maximum valuation of phase contrast in Kou, the method solving middle use iterative, constraints s.t ensures in estimating window, t Solution during secondary and the t-1 time iterative twines phase difference valuation and meets the requirement that difference is minimum, and K represents and estimates the big of window Little, it is ensured that (i, j) place's difference is minimum, and direction and the big I of K window carry out adaptive according to following rule for center phase point Should adjust: (i+n, j), (i, j+n), (i+n, j+n) three directions, the value of n ensures to search in the search the most respectively when difference becomes big Rope window is the window of a k × k；If changing the direction of search can not meet condition, then change the size of window: Round in expression；K=2 × K when difference diminishes, window expansion is twice and continues search for by expression.

Step 8, phase unwrapping model is changed into integer multiple difference estimation problem.

In conjunction with the characteristic of InSAR discrete phase matrix, in phase unwrapping, if phase contrast before and after only considering phase unwrapping It is two reciprocal processes that real part then phase unwrapping is wound around with phase place:

$\left\{\begin{array}{c}\mathrm{\&phi;}=\mathrm{\&psi;}+2\mathrm{\&pi;}(-k)\\ \mathrm{\&psi;}=\mathrm{\&phi;}+2\mathrm{\&pi;}k\end{array}\right.$

Wherein, φ is wound around phase place, and ψ is that solution twines rear absolute phase, and k is an integral multiple numerical value, then phase place (i, j) It is wound around phase contrast and absolute phase difference is defined as:

$\left\{\begin{array}{cc}{\mathrm{\&Delta;\&phi;}}_{i,j}^y=({\mathrm{\&psi;}}_{i,j+1}-{\mathrm{\&psi;}}_{i,j})+2\mathrm{\&pi;}(k_{i,j+1}-k_{i,j})& (i,j)\&Element;S_1\\ {\mathrm{\&Delta;\&phi;}}_{i,j}^x=({\mathrm{\&psi;}}_{i+1,j}-{\mathrm{\&psi;}}_{i,j})+2\mathrm{\&pi;}(k_{i+1,j}-k_{i,j})& (i,j)\&Element;S_2\end{array}\right.$

Image for a width M × N:

(i,j)∈S₁I=0,1 ..., M-2；J=0,1 ..., N-1

(i,j)∈S₂I=0,1 ..., M-1；J=0,1 ..., N-2

Absolute phase difference differs 2 π integral multiples with being wound around phase contrast:

$\left\{\begin{array}{cc}{\mathrm{\&Delta;\&phi;}}_{i,j}^y={\mathrm{\&Delta;\&psi;}}_{i,j}^y+2{\mathrm{\&pi;\&Delta;k}}_{i,j}^y& (i,j)\&Element;S_1\\ {\mathrm{\&Delta;\&phi;}}_{i,j}^x={\mathrm{\&Delta;\&psi;}}_{i,j}^x+2{\mathrm{\&pi;\&Delta;k}}_{i,j}^x& (i,j)\&Element;S_2\end{array}\right.$

Can be solved by above formula and be wound around phase contrast and solution and twine the discrete derivative deviation (after 2 π) of rear phase contrast, construction solution Twine longitudinal separation to orientation to the integer multiple of phase contrastWithFunctional equation be:

If it is rightWithCarry out estimating and integration can be obtained by the phase place after solution twines, asking of phase contrast will be solved The problem that topic is converted to solve integer multiple difference.

Utilize and precalculate the phase contrast statistical eigenfunction obtained, reconfigure solution and twine the function model of integer multiple difference Such as formula (13):

Wherein,WithFor utilize phase contrast statistical model estimate after recalculate the phase contrast integral multiple obtained Number, this valuation account for the situation of change of phase contrast, namely account for the change of orographic factor, and the solution obtained twines model more Accurately, reliably.

Step 9, utilize solution twine before and after's integer multiple difference minimum build solution twine model.

Least square phase restriction solution is twined model and replaces with the constrained minimization model of integer multiple difference, such as formula (14):

$\begin{array}{c}\min \left(J\right)=\{\underset{i=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left(P\right({\mathrm{\&Delta;\&phi;}}_x\left)\mathrm{\&Delta;}k\right)}_{i,j}^2+\underset{i=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left(P\right({\mathrm{\&Delta;\&phi;}}_y\left)\mathrm{\&Delta;}k\right)}_{i,j}^2\}\\ \begin{array}{ccc}s.t& \underset{K}{\&Sigma;}\&lsqb;\left(\right({\mathrm{Pk}}_{i,j+1}^y-{\mathrm{Pk}}_{i,j}^y)-({\mathrm{Pk}}_{i+1,j}^x-{\mathrm{Pk}}_{i,j}^x))({\mathrm{\&Delta;k}}^t-{\mathrm{\&Delta;k}}^{t-1})\&rsqb;<\mathrm{\&xi;}& t>1,(i,j)\&Element;S_{k\&times;k}\end{array}\end{array}---\left(14\right)$

Wherein, constraints s.t ensures estimating in window, and the t time and solution during the t-1 time iterative twine integral multiple Number difference meets the requirement that difference is minimum, and K represents the size estimating window；In discrete interference matrix,This solution twines model representation Littleization phase contrast statistical characteristics, i.e. twining phase contrast by solution and estimate the deviation being wound around phase contrast to be zero as far as possible.

The phase contrast field that discrete phase contrast departure is the absolute phase of irrotationality is obtained by formula (14), then along path Integration i.e. can obtain solution and twine phase place.Integral formula such as (15):

${\mathrm{\&psi;}}_{i,j}={\mathrm{\&phi;}}_0+2\mathrm{\&pi;}ent(\underset{i=0}{\overset{k-1}{\&Sigma;}}P\left({\mathrm{\&Delta;\&phi;}}_x\right){\mathrm{\&Delta;k}}_{i,j}^x+\underset{j=0}{\overset{k-1}{\&Sigma;}}P\left({\mathrm{\&Delta;\&phi;}}_y\right){\mathrm{\&Delta;k}}_{i,j}^y)---\left(15\right)$

Wherein, ent () represents rounding operation.

Step 10, solution to model calculate method.

The method utilizing discrete least-squares iteration solves solution and twines formula:

${\mathrm{Pk}}_{i,j}^{t+1}=({\mathrm{Pk}}_{i+1,j}^t+{\mathrm{Pk}}_{i-1,j}^t+{\mathrm{Pk}}_{i,j+1}^t+{\mathrm{Pk}}_{i,j-1}^t-{\mathrm{\&rho;}}_{i,j})/4$

Above formula represents that the t+1 time phase value is to be changed by the phase place of the t time to obtain, if (i, j) for estimating window center Pixel, ρ_i,jCan be obtained by following formula:

$\begin{array}{c}{\mathrm{\&rho;}}_{i,j}=k_{i,j}^{x}P\left({\mathrm{\&Delta;\&phi;}}_{i,j}^x\right)+k_{i,j}^{y}P\left({\mathrm{\&Delta;\&phi;}}_{i,j}^y\right)\\ =(k_{i+1,j}-2k_{i,j}+k_{i-1,j})P\left({\mathrm{\&Delta;\&phi;}}_{i,j}^x\right)+(k_{i,j+1}-2k_{i,j}+k_{i,j-1})P\left({\mathrm{\&Delta;\&phi;}}_{i,j}^y\right)\end{array}$

Phase restriction condition is all detected in iterative process each time, ifThen demonstrate,prove Bright for phase distortion, the maximum of probability phase contrast replacement in utilization estimation windowCarry out phase Position is replaced, it is ensured that the estimation window phase contrast statistical value of twice search is minimum, otherwise automatically changes window according to the method in step 7 The size and Orientation of mouth.

Step 11, when judging window search the most individual interferogram then disentanglement bundle, obtain disentanglement fruit.

By specific embodiment the present invention done above further describe it should be understood that, the most concrete Describing, should not be construed as the restriction to the spirit and scope of the invention, one of ordinary skilled in the art is reading this explanation The various amendments after book made above-described embodiment, broadly fall into the scope that the present invention is protected.

Claims (1)

1. an InSAR image phase unwrapping method based on phase contrast statistical model, it is characterised in that comprise the following steps:

Step one: assuming that pending interference image row, column pixel size is M × N, be divided into the partial estimation of k × k size Window, calculates phase place at the difference △ φ of ranks both direction, and to phase contrast in orientation to distance to derivation, obtain phase place Phase contrast derivative (phase frequency under discrete phase)

Wherein, △ x_i,j=x_i+1,j+x_i-1,j-2x_i,j, △ y_i,j=y_i,j+1+y_i,j-1-2y_i,jRepresent the interferogram phase along x, y direction respectively Potential difference value；It is the phase contrast average in x, y direction in k × k window respectively；With window center Phase contrast derivative standard deviation value；

Step 2: according to upper step phase contrast derivative value, calculates phase contrast statistical model in partial estimation window, and estimates local window The intra-orally distribution of shape high variable quantity；

Step 3: estimate Terrain Elevation variable quantity according to the relational expression of phase contrast frequency with terrain slope angle；

Step 4: obtain phase contrast statistical model in partial estimation window；

Step 5: utilize estimating window intra-orally shape variable quantity distribution function to smooth phase contrast in window, is estimating window Interior by distance to being multiplied by phase contrast statistical model obtains window the phase contrast taking topography variation amount into account to phase difference value with orientation Estimate model；

Step 6: construct InSAR image least square phase unwrapping model based on phase contrast statistical model；

Step 7: solution is twined longitudinal separation to and orientation to the integer multiple difference △ k of phase contrast_yWith △ k_xEstimate and amass Point, it is possible to obtain the phase contrast integer after solution twines, be converted to solve asking of integer multiple difference by the problem solving phase contrast Topic, the functional equation of structure integer multiple parameter is:WithIts Middle △ φ_yWith △ φ_xRepresent that phase point is wound around phase contrast,WithRepresent that solution twines rear phase contrast；

Step 8: least square phase unwrapping restricted model is replaced with the constrained minimization model of integer multiple difference, such as formula:

$\min \left(J\right)=\left\{\underset{i=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left(P\left({\mathrm{\&Delta;\&phi;}}_x\right)\mathrm{\&Delta;}k\right)}_{i,j}^2+\underset{i=0}{\overset{M-1}{\mathrm{\&Sigma;}}}\underset{j=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{\left(P\left({\mathrm{\&Delta;\&phi;}}_y\right)\mathrm{\&Delta;}k\right)}_{i,j}^2\right\}$

$\begin{array}{ccc}s.t& \underset{K}{\mathrm{\&Sigma;}}\&lsqb;\left(\right({\mathrm{Pk}}_{i,j+1}^y-{\mathrm{Pk}}_{i,j}^y)-({\mathrm{Pk}}_{i+1,j}^x-{\mathrm{Pk}}_{i,j}^x\left)\right)({\mathrm{\&Delta;k}}^t-{\mathrm{\&Delta;k}}^{t-1})\&rsqb;<\mathrm{\&xi;}& t>1,(i,j)\&Element;S_{k\&times;k}\end{array}$

Wherein, P (△ φ_x)、P(△φ_y) represent phase point row vector, column vector phase contrast statistical value, △ k^tWith △ k^t-1Represent The integer multiple of twice iteration is poor, and ξ is minimum；Constraints s.t represents in local window, the t time and the t-1 time iteration Solution when solving twines integer multiple difference and meets the requirement that difference is minimum, and K represents the size estimating window；

Step 9: utilize discrete least-squares iteration method to obtain solution and twine phase place.