CN101078769A - One-time all-polarization synthetic aperture radar image inverse method for digital elevation model - Google Patents

Abstract

The invention relates to a method of composing inversion digital elevation model of aperture radar image from single full-polarization by combining respective characteristics of radar polarization and radar angle-measurement technology. Polarization azimuth offset psi is estimated from polarization SAR data by using circular polarization method to build data relationship between polarization azimuth offset psi and range slope angle alpha and azimuth slope angle beta . Under the condition that there is only single full-polarization SAR image, suppose that imaging region satisfies conditions of surface covering homogeneity and little slope fluctuation and the technology from shadow to image in radar angle measurement is inducted to build data relationship of surface backward-scattering strength and surface geometric parameters (containing range slope angle alpha and azimuth slope angle beta ). In this way, by said data relationship range slope angle alpha and azimuth slope angle beta are calculated. Weight graph is extracted from cross-polarization response signal graph and discrete poisson equation about ground elevation is solved by weighted least square algorithm to obtain DEM in inversion.

Description

The method of one-time all-polarization synthetic aperture radar image inverse digital elevation model

Technical field

The invention belongs to space remote sensing and technical field of image processing, be specifically related to a kind of method of obtaining digital elevation model (DEM) by diameter radar image.

Background technology

(Synthetic Aperture Radar SAR) is a kind of airborne or satellite borne sensor of active to synthetic-aperture radar.Its ultimate principle is: satellite in orbital flight in diverse location periodically to same atural object launching electromagnetic wave pulse signal, receive echoed signal simultaneously.In some sense, can think to have extended the length of radar antenna, thereby improve resolution greatly.In addition, since synthetic-aperture radar also have round-the-clock, round-the-clock, be not subjected to advantages such as atmospheric propagation and climate effect, penetration power are strong, therefore civil and application military aspect is very extensive.

Digital elevation model (Digital Elevation Model), vehicle economy M.It is a kind of entity ground model of representing ground elevation with one group of orderly array of values form, is a digital terrain model branch of (Digital Terrain Model is called for short DTM).It is generally acknowledged, DTM describes the various geomorphologic factors that comprise elevation, as factors such as the gradient, aspect, slope change rates in the interior linearity and the space distribution of nonlinear combination, wherein DEM is the simple individual event numeral landform model of zeroth order, and other can derive from the basis of DEM as landforms characteristics such as the gradient, aspect and slope change rates.

At present, the method by SAR Image Acquisition digital elevation model (DEM) mainly contains following four kinds: radar angle measurement (Radarclinometry), radar stereopsis are measured (Radargrammetry), radar interference measures (SAP Interferometry) and polarization radar is measured (SAR Polarimetry).Wherein, the radar stereopsis is measured and radar interference is measured and must be used two width of cloth SAR images, and the former existence causes the problem of registration difficulty by different imagings visual angle, though and the latter's precision is high, the realization difficulty of correlation technique is big; Radar angle measurement only needs single width SAR image, direct backscatter intensity inverting DEM from the face of land, but constraint condition is too much in refutation process, and ignored the orientation to the gradient, so precision is not high; And the polarization radar measurement is a kind of the latest development direction in the D S AR imaging technique, this method is different from the three-dimensional imaging technology that radar interference is measured, not to utilize the phase differential of two width of cloth images to extract elevation information, but utilize the polarization integrated information of radar return to detect Terrain Elevation.

People such as Schuler studies show that, when the face of land of otherwise flat in the orientation when (being heading) goes up run-off the straight, the maximum position of its same polarization response can change, and promptly polarization orientation angle ψ will have a skew, and this phenomenon can be used for estimating terrain and rise and fall.Because the flight of single complete polarization can only obtain the orientation to the gradient by the polarization orientation angle estimation technique, yet also needs known distance to the gradient, could utilize the DEM of the whole scape image of least square technique construction of similar interference SAR phase unwrapping.Generally, need the quadrature flight of the reference DEM of survey region or another time complete polarization could obtain and the orientation to the distance of quadrature to the gradient.Because present polarization SAR data and with reference to the shortage of DEM, method above-mentioned often is not easy to realize.

Only obtaining the method for DEM in the prior art by the flight of single complete polarization, is the expression formula that obtains a relevant angle of range β and horizontal azimuth γ by the Eulerian angle conversion, obtains angle of range β in conjunction with the polarization orientation angle ψ that estimates again from polarization data.Again because the orientation to gradient SA with apart from representing with angle of range β and horizontal azimuth γ to gradient SR, obtains the DEM inverting with this.But the method need be used complicated digital image processing techniques when extracting horizontal azimuth γ.

Summary of the invention

The objective of the invention is to overcome the shortcoming of above-mentioned prior art, proposed a kind of inversion method that uses the digital elevation model (DEM) of single complete polarization in conjunction with polarization radar and radar angle measurement technology.

Method according to one-time all-polarization synthetic aperture radar image inverse digital elevation model of the present invention, comprise the steps: (1) use circular polarisation method, be offset from polarimetric synthetic aperture radar data estimation polarization orientation angle ψ, and go singular point to handle to ψ skew, set up polarization orientation angle ψ skew and well known to a person skilled in the art apart from the data relationship between angle of gradient β to angle of gradient α and orientation; (2) use in the radar angle measurement from the shade to the shape technology, set up face of land backscatter intensity with apart from the data relationship between angle of gradient β to angle of gradient α and orientation; (3) from above-mentioned steps (1) and (2) two data relation, computed range to gradient α and orientation to gradient β; (4) from polarimetric synthetic aperture radar data computation co-polarization response signal figure, from described co-polarization response signal figure, extract weight map; (5), find the solution the discrete Poisson equation of relevant ground elevation with the least-squares algorithm of weighting, thereby inverting obtains the digital elevation model of synthetic-aperture radar observation area from the result of above-mentioned steps (3) and (4).

The present invention be based on the polarization radar measuring method to the orientation to gradient sensitivity, and the radar angle measurement technology is adjusted the distance to the basis of gradient sensitivity, characteristics separately in conjunction with these two kinds of D S AR imaging techniques, from the method for single full polarimetric SAR inverting digital elevation model (DEM), thereby has easy realization, can satisfy the advantage of medium accuracy topographic mapping demand on a large scale.

Description of drawings

Fig. 1 detects synoptic diagram for singular point;

Fig. 2 is the geometric parameter model of bin;

Fig. 3 is a calculation procedure FB(flow block) of the present invention;

Fig. 4 is the SIR-C data (C-band Pauli decomposes the RGB image) in area, province Urumchi, Chinese Xinjiang;

Fig. 5 is that orientation, area, Urumchi is to slope map;

Fig. 6 is that area, Urumchi distance is to slope map;

Fig. 7 is the weight quality figure in area, Urumchi;

The DEM figure that Fig. 8 is finally inversed by by the present invention for the area, Urumchi;

Fig. 9 is the reference DEM figure in area, Urumchi;

Figure 10 is the SIR-C data (L-band Pauli decomposes the RGB image) in area, Russian look rib lattice river;

Figure 11 is that the orientation in area, Russian look rib lattice river is to the slope map picture;

Figure 12 is that the distance in area, Russian look rib lattice river is to the slope map picture;

Figure 13 is the weight quality figure in area, Russian look rib lattice river;

Figure 14 is the DEM figure of area, Russian look rib lattice river by inverting of the present invention;

Figure 15 is the contour map of the inverting DEM in area, Russian look rib lattice river.

Embodiment

With reference to the accompanying drawings the present invention is specifically described.The method according to this invention at first, utilizes the circular polarisation method to estimate polarization orientation angle ψ skew from the polarization SAR image, sets up polarization orientation angle ψ skew and terrain slope, promptly apart from the data relationship between angle of gradient β to angle of gradient α and orientation.Under the situation of having only the single full polarimetric SAR, suppose that imaging region satisfies homogeney and the little condition of gradient fluctuating that the face of land covers, introduce (the Shape from Shading) technology from the shade to the shape in the radar angle measurement, the data relationship between the backscatter intensity of setting up the face of land and the face of land geometric parameter (comprise apart to angle of gradient α and orientation to angle of gradient β).Like this, by above-mentioned two data relation, just can directly calculate apart to gradient α and orientation to gradient β.At last, from co-polarization response signal figure, extract weight map, find the solution the discrete Poisson equation of relevant ground elevation with the least-squares algorithm of weighting, thereby inverting obtains DEM.

Fig. 1 detects synoptic diagram for singular point; Fig. 2 is the geometric parameter model of bin; Fig. 3 is a calculation procedure FB(flow block) of the present invention.1,2 and 3 pair of content of the present invention specifically describes with reference to the accompanying drawings.

I. the estimation and the optimization of polarization orientation angle ψ skew

(1) input is used to generate polarization matrix [C3] or [T3] data of DEM, and radar visual angle φ and image resolution ratio (comprise that distance is to resolution R _dWith the orientation to resolution R _a).

(2) polarization matrix [C3] or [T3] to input carries out polarization SAR filtering, can select Boxcar and Refined Lee polarization filtering;

(3) obtain filtered polarization matrix [C3] or [T3] by above step, utilize the polarization orientation angle method of estimation to carry out preliminary polarization orientation angle and estimate.Krogager and Czyz propose to estimate the position angle with dextrorotation and left-hand circular polarization.This method is further revised and is refined by people such as Lee and obtains the circular polarisation method, and it is as follows specifically to derive:

Polarization data generally provides two kinds of polarization type: multiple (SLC) data of haplopia and look multiple (MLC) data more.Look complex data owing to use the haplopia complex data to derive, so we are that example is carried out preliminary polarization orientation angle estimation with the haplopia complex data at this more.

The polarization SAR data can be used scattering matrix (1) expression

$S=\left[\begin{array}{cc}{S}_{\mathrm{hh}}& S_{\mathrm{hv}}\\ S_{\mathrm{vh}}& S_{\mathrm{vv}}\end{array}\right]---\left(1\right)$

Change emission or reception.

Because the SAR satellite all is single based system basically at present, promptly transmitting and receiving antenna is the same antenna, satisfies the antenna reciprocity condition, i.e. S _Hv=S _VhAzimuthal rotation obtains by following formula

$\tilde{S}=\left[\begin{array}{cc}\cos \left(\mathrm{\ψ}\right)& \sin \left(\mathrm{\ψ}\right)\\ -\sin \left(\mathrm{\ψ}\right)& \cos \left(\mathrm{\ψ}\right)\end{array}\right]\left[\begin{array}{cc}{S}_{\mathrm{hh}}& S_{\mathrm{hv}}\\ S_{\mathrm{vh}}& S_{\mathrm{vv}}\end{array}\right]\left[\begin{array}{cc}\cos \left(\mathrm{\ψ}\right)& -\sin \left(\mathrm{\ψ}\right)\\ \sin \left(\mathrm{\ψ}\right)& \cos \left(\mathrm{\ψ}\right)\end{array}\right]---\left(2\right)$

In expression formula (2) formula above the matrix S～represent postrotational matrix.

Each circular polarisation component can obtain from scattering matrix:

S _rr＝(S _hh-S _vv+i2S _hv)/2

S _ll＝(S _vv-S _hh+i2S _hv)/2 (3)

S _rl＝i(S _hh+S _vv)/2

Azimuthal rotation can be used expression formula (2) formula and be obtained, and after the arrangement is

${\tilde{S}}_{\mathrm{rr}}=S_{\mathrm{rr}}{e}^{-i2\mathrm{\ψ}}$

${\tilde{S}}_{\mathrm{ll}}=S_{\mathrm{ll}}{e}^{i2\mathrm{\ψ}}---\left(4\right)$

${\tilde{S}}_{\mathrm{rl}}=S_{\mathrm{rl}}$

Define a circular polarisation vector

$c=\left[\begin{array}{c}{S}_{\mathrm{rr}}\\ \sqrt{2}{S}_{\mathrm{rl}}\\ S_{\mathrm{ll}}\end{array}\right]---\left(5\right)$

The covariance matrix G of circular polarisation base can obtain from vectorial c

$G=\⟨{\mathrm{cc}}^{*T}\⟩=\left[\begin{array}{ccc}\⟨{\left|S_{\mathrm{rr}}\right|}^2\⟩& \sqrt{2}\⟨\left(S_{\mathrm{rr}}{S}_{\mathrm{rl}}^{*}\right)\⟩& \⟨\left(S_{\mathrm{rr}}{S}_{\mathrm{ll}}^{*}\right)\⟩\\ \sqrt{2}\⟨\left(S_{\mathrm{rl}}{S}_{\mathrm{rr}}^{*}\right)\⟩& 2\⟨{\left|S_{\mathrm{rl}}\right|}^2\⟩& \sqrt{2}\⟨\left(S_{\mathrm{rl}}{S}_{\mathrm{ll}}^{*}\right)\⟩\\ \⟨\left(S_{\mathrm{ll}}{S}_{\mathrm{rr}}^{*}\right)\⟩& \sqrt{2}\⟨\left(S_{\mathrm{ll}}{S}_{\mathrm{rl}}^{*}\right)\⟩& \⟨{\left|S_{\mathrm{ll}}\right|}^2\⟩\end{array}\right]---\left(6\right)$

Use expression formula (4), the circular polarisation covariance matrix that rotates through For

$\tilde{G}=\left[\begin{array}{ccc}\&lang;{\left|S_{\mathrm{rr}}\right|}^2\&rang;& \sqrt{2}\&lang;\left(S_{\mathrm{rr}}{S}_{\mathrm{rl}}^{*}\right)e^{-\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">i</figure-callout>}2\mathrm{\&psi;}}\&rang;& \&lang;\left(S_{\mathrm{rr}}{S}_{\mathrm{ll}}^{*}\right)e^{-\mathrm{<figure-callout\; id="4"\; label="i"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">i</figure-callout>}4\mathrm{\&psi;}}\&rang;\\ \sqrt{2}\&lang;\left(S_{\mathrm{rl}}{S}_{\mathrm{rr}}^{*}\right)e^{i2\mathrm{\&psi;}}\&rang;& 2\&lang;{\left|S_{\mathrm{rl}}\right|}^2\&rang;& \sqrt{2}\&lang;\left(S_{\mathrm{rl}}{S}_{\mathrm{ll}}^{*}\right)e^{-\mathrm{<figure-callout\; id="2"\; label="i"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">i</figure-callout>}2\mathrm{\&psi;}}\&rang;\\ \&lang;\left(S_{\mathrm{ll}}{S}_{\mathrm{rr}}^{*}\right)e^{i4\mathrm{\&psi;}}\&rang;& \sqrt{2}\&lang;\left(S_{\mathrm{ll}}{S}_{\mathrm{rl}}^{*}\right)e^{i2\mathrm{\&psi;}}\&rang;& \&lang;{\left|S_{\mathrm{ll}}\right|}^2\&rang;\end{array}\right]---\left(7\right)$

Obviously, for the circular polarisation covariance matrix

, azimuthal variation only influences the phase place of off diagonal element item.The right side-the right side and a left side-left circular polarisation item, promptly (1,3) item can be used to estimate the position angle.If azimuthal skew mainly causes to the gradient by the orientation, so

$\&lang;{\tilde{S}}_{\mathrm{rr}}{\tilde{S}}_{\mathrm{ll}}^{*}\&rang;=\&lang;S_{\mathrm{rr}}{S}_{\mathrm{ll}}^{*}\&rang;e^{-\mathrm{<figure-callout\; id="4"\; label="i"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">i</figure-callout>}4\mathrm{\&psi;}}---\left(8\right)$

By following formula as can be known, have only as＜S _RrS _Ll ^*When being real number, just can not destroy and phase term e ^{-i4 ψ}Relevant position angle.Because the most of atural objects of nature have reflection symmetry, and can prove, for reflective symmetry medium,＜S _RrS _Ll ^*Be real number, intersection and co-polarization continuous item are 0.With formula (3) substitution＜S _RrS _Ll ^*, and establish and comprise S _HvItem (except＜S _Hv| ²) be 0, then have

$\&lang;S_{\mathrm{rr}}{S}_{\mathrm{ll}}^{*}\&rang;=(-\&lang;{|S_{\mathrm{hh}}-S_{\mathrm{vv}}|}^2\&rang;+4\&lang;{\left|S_{\mathrm{hv}}\right|}^2\&rang;)/4---\left(9\right)$

This is a real number, therefore＜and S _RrS _Ll ^*Argument be 0 or π.So, S _HhAnd S _VvBetween phase differential can not produce deviation to Azimuth Estimation.

The scope that the factor 4 ψ of expression formula (8) have limited the ψ angle is [π/4, π/4].In order to derive general expression formula, we use linear basis representation

$\&lang;{\tilde{S}}_{\mathrm{rr}}{\tilde{S}}_{\mathrm{ll}}^{*}\&rang;=\frac{1}{4}\{\&lang;-{|{\tilde{S}}_{\mathrm{hh}}-{\tilde{S}}_{\mathrm{vv}}|}^2+4{\left|{\tilde{S}}_{\mathrm{hv}}\right|}^2\&rang;-\mathrm{<figure-callout\; id="4"\; label="i"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">i</figure-callout>}4\mathrm{Re}(\&lang;({\tilde{S}}_{\mathrm{hh}}-{\tilde{S}}_{\mathrm{vv}}){\tilde{S}}_{\mathrm{hv}}^{*}\&rang;)\}---\left(10\right)$

To cause deviation if directly use following formula, because for the reflective symmetry medium,

Usually greater than

Denominator is a negative.Therefore, when molecule near 0 the time, arc-tangent value is approaching ± π, the position angle will be ± π/4, rather than approach 0 as it is right.For position angle and corresponding orientation are complementary to the gradient, deviation must remove by adding π.Finally, can obtain circular polarisation is estimated as

$\mathrm{\&psi;}=\left\{\begin{array}{cc}\mathrm{\&eta;}& \mathrm{\&eta;}\&le;\mathrm{\&pi;}/4\\ \mathrm{\&eta;}-\mathrm{\&pi;}/2& \mathrm{\&eta;}>\mathrm{\&pi;}/4\end{array}\right.---\left(11\right)$

Wherein

$\mathrm{\&eta;}=\frac{1}{4}[{\tan }^{-1}\left(\frac{-4\mathrm{Re}(\&lang;({\tilde{S}}_{\mathrm{hh}}-{\tilde{S}}_{\mathrm{vv}}){\tilde{S}}_{\mathrm{hv}}\&rang;)}{-\&lang;{|{\tilde{S}}_{\mathrm{hh}}-{\tilde{S}}_{\mathrm{vv}}|}^2\&rang;+4\&lang;{\left|{\tilde{S}}_{\mathrm{hv}}\right|}^2\&rang;}\right)+\mathrm{\&pi;}]---\left(12\right)$

The computer capacity of arc tangent is at [π, π] in the expression formula (12).

(4) estimated accuracy of polarization orientation angle depends on whether the position angle skew of main scatterer is only caused to the gradient by the orientation.The mixing of the scatterer of different scattering mechanisms is equivalent to the stack of scattering (or relevant) matrix, when the identical position angle of rotation, can not bring error, yet, if have the scattering of some to have because the position angle translation that branch or other factors cause so just may have problems.At this moment often because the depolarization effect is obvious or peak point is not unique, the polarization orientation angle ψ skew that therefore estimates may the angle of arrival singular point.Especially, the present invention has adopted and has removed the singular point treatment technology.Its concrete steps are as follows:

A) an angle change threshold T is set;

B) singularity of detection pixel.The order of pressing Row Column scans each pixel successively in entire image.As shown in Figure 1,0 expression current pixel, 1,2,3,4 expressions are adjacent four pixels, if when the angle difference of arbitrary pixel is greater than given threshold value T in the current pixel 0 and 1,4, judge that then current pixel point is a singular point.Do like this be in order to ensure this singular point in singular regions inside, rather than the non-singular point at singular regions edge.

C) if current pixel point is a non-singular point, then return step (b); Otherwise, enter step (d).

D) intermediate value with all pixel values in the given window size (as 3*3) replaces current unusual pixel.

E) repeat (b)-(d) step 3 and time get final product, obtain the polarization orientation angle behind the singular point at last.

II. the foundation of image/face of land relation

(1) polarization matrix [C3] or [T3] the extracting data HH polarization intensity figure that imports from step I (1), and carry out SAR image filtering (as Lee, Frost, Gamma etc.).But note not using enhancement filter, because its can retention point target, and the existing of point target (as isolated trees etc.) can't truly be reflected the true altitude on the face of land, is disadvantageous for three-dimensional elevation inverting.

(2) mean value of view picture HH polarization intensity figure after the calculation of filtered is as level land backscatter intensity I ₀Estimated value;

(3) Lee and Pottier independently derive polarization orientation angle ψ and terrain slope respectively and the radar visual angle has following relation

$\tan \mathrm{\&psi;}=\frac{\tan \mathrm{\&beta;}}{-\tan \mathrm{\&alpha;}\cos \mathrm{\&phi;}+\sin \mathrm{\&phi;}}---\left(13\right)$

Wherein, φ is radar visual angle (input in step I (1)), and tan β is that the orientation is to gradient S _A(subscript A represent orientation to), tan α are that distance is to gradient S _R(subscript R represent distance to).Because two unknown quantitys are arranged, and (α β) with an equation (8), obtain grade information, needs another relational expression.Especially, the present invention adopts (Shapefrom Shading) technology from the shade to the shape in the radar angle measurement method to obtain another about the expression formula between image/face of land.

Being used for optical imagery from the shade to the shape at first, is a kind of technology of recovering the respective objects shape under the known situation of reflective function from the reflection strength image of working as.R.L.Widley at first will be applied in the radar image from the shade to the shape technology, and rename and be radar angle measurement.Its principle is to estimate corresponding face of land terrain parameter from the backscatter intensity of SAR image, uses single width SAR image to disclose surface slope information fully.

From the shade to the shape technology, make us can estimate the geometric parameter on the corresponding face of land, need assume picture overlay area, the face of land homogeneity and the gradient and rise and fall mild by the backscatter intensity of atural object.For the geometric parameter on the face of land, we introduce the notion of bin (Parcel), and it is defined as the face of land that a radar resolution element is intercepted, and are imaged in the pixel.

Fig. 2 is radar image sampled point ground (when the gradient is not very big) corresponding face of land shape approximation figure when comparatively smooth.

At this moment, though L (α) is vertical with the unit area and the oblique distance of getting along well that L (β) constitutes, since in easy flights, can be similar to and think that it is exactly an area A.For each image-generating unit, comprise that three how much imaging parameters: the angle of depression θ of radar incident electromagnetic wave (φ is surplus mutually with the radar visual angle) distances are to resolution R _d, the orientation is to resolution R _a, therefore,

L(α)＝R _d sin(θ)/sin|θ-α| (14)

L(β)＝R _a/cos(β) (15)

A(α，β)＝L(α)L(β) (16)

Wherein φ, R _d, R _aInput in above-mentioned steps I (1).In addition, each pixel intensity not only depends on the gradient on imaging ground, and depends on its face of land cover type.Reflection model is arranged in optical system, backscattering model is then arranged in radar system.In order to set up the relation between ground geometric parameter and the SAR image intensity, we adopt simple and effective experience backscattering model--lambert's backscattering model.Can set up relation between image/ground by lambert hypothesis:

I _r(α，β)＝Kσ ₀R _dR _acos ²(θ-α)cos ²βA(α，β) (17)

Wherein, I _rBe backscatter intensity, α be the distance to the angle of gradient, β be the orientation to the angle of gradient, θ is the angle of depression of radar incident electromagnetic wave, K is the radar calibration constant, σ ₀Be the backscattering coefficient of mulching material, A is the floor area of a pixel correspondence on the SAR image.This relational expression also is the basis of radar angle measurement.

Because unknown quantity (K, the σ of each pixel ₀, α, β) numerous, in order to solve terrain slope from expression formula (10), we suppose the backscatter intensity I for a smooth ground ₀, floor type is identical, so σ ₀Be constant in quite on a large scale, calculate the backscatter intensity I on level land so ₀The time α and β be 0, promptly

I ₀＝I _r(0，0)＝Kσ ₀R _dR _acos ²(θ) (18)

The level land backscatter intensity I of this moment ₀Can think constant.Under normal conditions, we can be used as I to the assembly average of a width of cloth image ₀, based on this, the proportionate relationship below setting up, be designated as R (α, β):

$R(\mathrm{\&alpha;},\mathrm{\&beta;})=\frac{I_r(\mathrm{\&alpha;},\mathrm{\&beta;})}{I_0}=\frac{\sin \mathrm{\&theta;}{\cos }^2(\mathrm{\&theta;}-\mathrm{\&alpha;})\cos \mathrm{\&beta;}}{{\cos }^2\left(\mathrm{\&theta;}\right)\sin \left(\right|\mathrm{\&theta;}-\mathrm{\&alpha;}\left|\right)}---\left(19\right)$

By finding the solution following formula, we have obtained the expression formula of distance to gradient α:

$\mathrm{\&alpha;}=\mathrm{\&theta;}-\arcsin ((\sqrt{R^{\&prime;2}+4}-R^{\&prime;})/2)---\left(20\right)$

Wherein,

$R^{\&prime;}=\frac{R(\mathrm{\&alpha;},\mathrm{\&beta;}){\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">cos</figure-callout>}}^2\mathrm{\&theta;}}{\sin \mathrm{\&theta;}\cos \mathrm{\&beta;}}---\left(21\right)$

Therefore, in conjunction with the characteristics separately of polarization radar and radar angle measurement, from top two data relations (13) and (20), (α β), thereby obtains the orientation to gradient S just can directly to calculate the surface slope angle _AWith the distance to gradient S _R

III. by surface slope S (x, y) inverting DEM

(1) being similar to phase unwrapping in the interference SAR by terrain slope inverting DEM, is the separating of discrete Poisson (Possion) equation of seeking to have the Newman edge-restraint condition on M * N rectangular node zone, and this Poisson equation is expressed as

 ²h(x，y)＝ρ(x，y) (22)

Wherein, (x y) is ground elevation,  to h ²Be Laplace operator  ²/  x ²+  ²/  y ²Use h _IjBe illustrated in (i, the height value of j) locating, source function ρ _IjThen be by input gradient data S _IjThe surface curvature value that obtains,

${\mathrm{\&rho;}}_{\mathrm{ij}}^R=S_{\mathrm{ij}}^R-S_{i-1,j}^R,{\mathrm{\&rho;}}_{\mathrm{ij}}^A=S_{\mathrm{ij}}^A-S_{i,j-1}^A(0\&le;i\&le;M,0\&le;j\&le;N)---\left(23\right)$

In order to obtain ground elevation, must minimize h _IjDifference between the local derviation is used following sum formula

${\mathrm{\&chi;}}^2=\underset{i,j}{\mathrm{\&Sigma;}}{(h_{\mathrm{ij}}-h_{i-1,j}-S_{\mathrm{ij}}^R)}^2+\underset{i,j}{\mathrm{\&Sigma;}}{(h_{\mathrm{ij}}-h_{i,j-1}-S_{\mathrm{ij}}^A)}^2---\left(24\right)$

To h _IjDifferentiate, and make its result equal 0, then obtain

(h _i-1，j-2h _ij+h _i-1，j)+(h _i，j+1-2h _ij+h _i，j+1)＝ρ _ij (25)

The Poisson equation that Here it is will find the solution.Wherein, Li San source function ρ _Ij

${\mathrm{\&rho;}}_{\mathrm{ij}}=(S_{i+1,j}^R-S_{\mathrm{ij}}^R)+(S_{i,j+1}^A-S_{\mathrm{ij}}^A)---\left(26\right)$

The method for solving of prior art has alternating direction implicit iteration (ADI) algorithm and the complete multi grid of no weight (FMG) algorithm.

FMG is a kind of fast method that solves partial differential equation (PDE) on than the macrolattice basis of classics, carries out the calculating of Gauss-Seidel relaxation method on thicker, littler grid.Usually, for solving discrete PDE problem on N * N grid, FMG is the same fast with the direct method speed based on Fourier transform, and to the size of grid without limits.Its key is and can be converted to the high frequency composition to the low-frequency component of error, and radio-frequency component can be removed with the Gauss-Seidel relaxation method soon like this, and this is to realize by this problem is transferred in the thicker grid.

Though there is not the least square method fast operation of weight, it has an obvious defects, i.e. be the inconsistent district of passing the gradient during its inverting, rather than walk around, so usually can cause the error of topographical surface.This defective can be remedied by introducing weight.Especially, the present invention adopts the FMG algorithm of weighting.In general, weight data w _{I, j}Span be [0,1] because polarization SAR does not exist phase dry mass figure in the interference SAR phase unwrapping as weight map, so we define similar weight quality figure from co-polarization signal response figure

$w_{i,j}=\frac{{\mathrm{\&sigma;}}_{\max }(0,{\mathrm{\&psi;}}^{\&prime;})-{\mathrm{\&sigma;}}_{\min }(0,{\mathrm{\&psi;}}^{\&prime;})}{{\mathrm{\&sigma;}}_{\max }(0,{\mathrm{\&psi;}}^{\&prime;})}---\left(27\right)$

Wherein, σ _Max(0, ψ ') and σ _Min(0, ψ ') expression considered pixel (i, j) maximum/minimum response of the co-polarization signal under linear polarization base (χ=0) situation.

The backscattering cross of target (or received power) is the function of incident wave polarized state, promptly is the polarization ellipse rate of incident wave and the function of polarization orientation angle.Polarization response is called the co-polarization response when the polarized state that transmits and receives antenna is consistent, and the polarization response when their polarized state quadratures is called the cross polarization response.Here, only consider the co-polarization response, the backscattering cross of target can be expressed as so

$\mathrm{\&sigma;}(\mathrm{\&chi;},\mathrm{\&psi;})={\left[\begin{array}{c}1\\ \mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">cos</figure-callout>}2\mathrm{\&chi;}\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">cos</figure-callout>}2{\mathrm{\&psi;}}^{\&prime;}\\ \mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">cos</figure-callout>}2\mathrm{\&chi;}\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">sin</figure-callout>}2{\mathrm{\&psi;}}^{\&prime;}\\ \mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">sin</figure-callout>}2\mathrm{\&chi;}\end{array}\right]}^{\mathrm{<figure-callout\; id="1"\; label="T\; K"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">T</figure-callout>}}K\left[\begin{array}{c}1\\ \mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">cos</figure-callout>}2\mathrm{\&chi;}\mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">cos</figure-callout>}2{\mathrm{\&psi;}}^{\&prime;}\\ \mathrm{<figure-callout\; id="2"\; label="cos"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">cos</figure-callout>}2\mathrm{\&chi;}\mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">sin</figure-callout>}2{\mathrm{\&psi;}}^{\&prime;}\\ \mathrm{<figure-callout\; id="2"\; label="sin"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">sin</figure-callout>}2\mathrm{\&chi;}\end{array}\right]---\left(28\right)$

Wherein, K is the Stokes matrix; χ is the ellipticity angle, and span exists

ψ ' is a polarization orientation angle, and span exists

Formula (28) travels through whole (χ, ψ ') plane, can obtain the co-polarization response diagram of target.When target scattering during near single surface scattering, w _{I, j}Big more; When target scattering near volume scattering, in other words, when depolarisation effect strengthens, can cause increasing of co-polarization signal response figure substrate, so w _{I, j}More little.

Like this, having no right least square problem just is converted to the weighted least-squares problem.Consider following weighted least-squares problem:

$\min \{\underset{i,j}{\mathrm{\&Sigma;}}{w}_{\mathrm{ij}}^R{(h_{i+1,j}-h_{i,j}-S_{i,j}^R)}^2+\underset{i,j}{\mathrm{\&Sigma;}}{w}_{\mathrm{ij}}^A{(h_{i,j+1}-h_{i,j}-S_{i,j}^A)}^2\}---\left(29\right)$

Gradient weight w _Ij ^RAnd w _Ij ^ABe defined as follows:

$\begin{array}{cc}{w}_{\mathrm{ij}}^R=\min (w_{i+1,\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="A20061008127700251.png,A20061008127700261.png"\; state="\{\{state\}\}">j</figure-callout>}}^2,w_{i,j}^2)& w_{\mathrm{ij}}^A=\min (w_{i,j+1}^2,w_{i,j}^2)\end{array}---\left(30\right)$

So, the least square solution h of this problem _{I, j}Be defined as

${\widehat{\mathrm{\&rho;}}}_{\mathrm{ij}}=w_{i,j}^R(h_{i+1,j}-h_{i,j})-w_{i-1,j}^R(h_{i,j}-h_{i-1,j})$

$+w_{\mathrm{ij}}^A(h_{i,j+1}-h_{i,j})-w_{i,j-1}^A(h_{i,j}-h_{i,j-1})---\left(31\right)$

Wherein,

Be weighted phases Laplce operand,

${\widehat{\mathrm{\&rho;}}}_{\mathrm{ij}}=w_{\mathrm{ij}}^R{S}_{\mathrm{ij}}^R-w_{i-1,j}^R{S}_{i-1,j}^R+w_{\mathrm{ij}}^A{S}_{i,j}^A-w_{i,j-1}^A{S}_{i,j-1}^A---\left(32\right)$

Classic Gauss-Seidel relaxation method is passed through equation

$h_{i,j}=\frac{w_{\mathrm{ij}}^R{h}_{i+1,j}+w_{i-1,j}^R{h}_{i-1,j}+w_{\mathrm{ij}}^A{h}_{i,j+1}+w_{i,j-1}^A{h}_{i,j-1}-{\widehat{\mathrm{\&rho;}}}_{\mathrm{ij}}}{v_{i,j}}---\left(33\right)$

Iterative h repeatedly _{I, j}

In the expression formula (33), $v_{i,j}=w_{\mathrm{ij}}^R+w_{i-1,j}^R+w_{\mathrm{ij}}^A+w_{i,j-1}^A---\left(34\right)$

(2) use distance that step I (1) obtains to resolution R _dWith the orientation to resolution R _a, the orientation that obtains of Step II (3) is to gradient S _AWith the distance to gradient S _RWith the weight quality figure W that Step II I (1) obtains, use weighting FMG algorithm just can obtain DEM.

The elevation of weighting FMG algorithm by gradient data fitting topographical surface distinguished by the mode of 0 weight is set data " damages ", prevents error transfer, utilizes the Gauss-Seidel on the grid of multiresolution to relax accelerating convergence then.

Experiment 1: area, Urumchi, Xinjiang

The C-band data that on October in 1994 8, SIR-C observed Chinese Tianshan Area that adopt California Inst Tech USA jet propulsion laboratory (JPL) to provide, concrete parameter is as follows:

Experiment place: Chinese Tianshan Mountains

Imaging time: 1994/10/08 22:06:45.647

Picture centre latitude (degree): 43.6740608

Picture centre longitude (degree): 88.1214600

Picture centre incident angle (degree): 24.098

Radar frequency (GHz): 5.298

Product type: haplopia multiple (SLC)

Line-spacing (rice): 5.0520549

Pixel separation (rice): 6.6624818

The pixel count of every capable image: 1100

Picturedeep: 9896

Map projection's descriptor: oblique distance

Orbital direction: rail lift

Fig. 4 is that the size that intercepts from former figure is the Pauli decomposition composite diagram (be that R is HH-VV, G is HV, and B is HH+VV) of 273 * 367 pixels.Distance is 16.3177m to resolution, and the orientation is 15.1562m (processing is looked to having done three in the orientation) to resolution, and its center is that (88.0322E 43.5610N), is positioned at area, province Urumchi, Chinese Xinjiang.

From Fig. 4, Fig. 8 and Fig. 9 as can be known, the DEM of Fig. 8 inverting and the main atural object among Fig. 4 move towards basically identical, this also is confirmed in the reference DEM of Fig. 8 once more, and the light and shade texture of four zanjons among Fig. 8 among four of the left side zanjons and Fig. 9 and trend also meet very much.

In order to estimate the DEM inversion accuracy more objectively, the DEM of inverting among Fig. 8 is registrated on the reference DEM among Fig. 9, evenly choose 20 DEM reference mark in both common region, the detailed precision analysis that obtains is as shown in table 1:

The precision analysis at reference mark among the table 1 inverting DEM

Index	Latitude (N)	Longitude (E)	With reference to height value (m)	Measure height value (m)	Elevation deviation (m)
						01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20	43:30:56 43:30:35 43:31:40 43:30:02 43:30:37 43:31:10 43:30:54 43:30:41 43:30:26 43:31:14 43:30:59 43:31:03 43:30:02 43:30:34 43:31:02 43:30:04 43:30:20 43:31:03 43:31:31 43:30:33	88:02:04 88:02:31 88:03:01 88:02:18 88:03:27 88:02:44 88:02:48 88:02:53 88:01:50 88:01:35 88:03:17 88:01:48 88:02:55 88:02:36 88:01:45 88:02:36 88:03:26 88:01:29 88:02:18 88:02:24	1220.60 1220.60 1264.90 1190.10 1240.00 1249.70 1234.90 1230.00 1199.50 1206.40 1248.70 1215.20 1217.50 1224.00 1216.10 1200.40 1236.80 1207.30 1239.90 1217.50	1231.15 1229.43 1227.83 1217.19 1215.17 1241.82 1229.86 1217.21 1229.58 1232.83 1221.85 1232.85 1208.94 1225.64 1234.22 1211.91 1214.10 1232.99 1225.36 1229.46	-10.55 -8.83 37.07 -27.09 24.83 7.88 5.04 12.79 -30.08 -26.43 26.85 -17.65 8.56 -1.64 -18.12 -11.51 22.70 -25.69 14.54 -11.96

Average :-0.91 root-mean-square error: 20.29

Calculate the dispersed elevation deviation according to the reference mark of 20 of even distributions and be-0.91m, root-mean-square error is 20.29m.This precision can satisfy the demand of topographic mapping on a large scale fully.By table 1 as seen, the measured elevation that obtains of the present invention is compared with the reference elevation and is existed certain error, one of its source of error is from measuring DEM and with reference to the registration error between the DEM, the 2nd, because polarization data is gathered in C-band, its short wavelength is limited to vegetation coverture penetration depth, the elevation on the feasible elevation of measuring and the actual face of land is deviation to some extent, uses the longer L-band of wavelength then to be expected to improve the DEM precision.

Experiment 2: area, Russian look rib lattice river

The L-band data that on October in 1994 7, SIR-C observed area, Russian look rib lattice river that adopt California Inst Tech USA jet propulsion laboratory (JPL) to provide, concrete parameter is as follows:

Experiment place: area, Russian look rib lattice river

Imaging time: 1994/10/07 22:33:30.642

Picture centre latitude (degree): 52.1651077

Picture centre longitude (degree): 106.6797943

Picture centre incident angle (degree): 24.569

Radar frequency (GHz): 1.249

Product type: haplopia multiple (SLC)

Line-spacing (rice): 5.0520549

Pixel separation (rice): 6.6624818

The pixel count of every capable image: 1100

Picturedeep: 9030

Map projection's descriptor: oblique distance

Orbital direction: rail lift

Figure 10 is that the size that intercepts from former figure is the Pauli decomposition composite diagram (be that R is HH-VV, G is HV, and B is HH+VV) of 562 * 380 pixels.Distance is 16.3177m to resolution, and the orientation is 15.1562m (processing is looked to having done three in the orientation) to resolution, and its center is that (107.0003E 52.2443N), is positioned at area, Russian look rib lattice river.

Because this area can't obtain the higher reference DEM that can compare with inverting DEM of precision in the Russia area.But, can find out easily very intuitively that from Figure 10 and Figure 12 DEM that Figure 14 inverting obtains and image 10 mesorelief rise and fall and coincide finely, some trickleer topographic relieies of Figure 10 lower right-most portion change, and also in inverting DEM Figure 14 good embodiment are arranged.

Claims (5)

1, a kind of method of one-time all-polarization synthetic aperture radar image inverse digital elevation model, this method utilized polarization orientation angle ψ skew, distance to angle of gradient α and orientation to angle of gradient β, it is characterized in that, comprise the steps:

(1) use the circular polarisation method, from polarimetric synthetic aperture radar data estimation polarization orientation angle ψ skew, and to ψ skew goes singular point to handle, set up polarization orientation angle ψ skew with apart from the data relationship between angle of gradient β to angle of gradient α and orientation;

(2) use in the radar angle measurement from the shade to the shape technology, set up face of land backscatter intensity with apart from the data relationship between angle of gradient β to angle of gradient α and orientation;

(3) from above-mentioned steps (1) and (2) two data relation, computed range to gradient α and orientation to gradient β;

(4) from polarimetric synthetic aperture radar data computation co-polarization response signal figure, from described co-polarization response signal figure, extract weight map;

(5), find the solution the discrete Poisson equation of relevant ground elevation with the least-squares algorithm of weighting, thereby inverting obtains the digital elevation model of synthetic-aperture radar observation area from the result of above-mentioned steps (3) and (4).

2, method according to claim 1 is characterized in that, goes singular point to handle to the estimation of polarization orientation angle ψ skew, and wherein, what polarization orientation angle ψ was offset is estimated as

$\mathrm{\&psi;}=\left\{\begin{array}{cc}\mathrm{\&eta;}& \mathrm{\&eta;}\&le;\mathrm{\&pi;}/4\\ \mathrm{\&eta;}-\mathrm{\&pi;}/2& \mathrm{\&eta;}>\mathrm{\&pi;}/4\end{array}\right.---\left(11\right)$

Wherein

$\mathrm{\&eta;}=\frac{1}{4}[{\tan }^{-1}\left(\frac{-4\mathrm{Re}(\&lang;({\tilde{S}}_{\mathrm{hh}}-{\tilde{S}}_{\mathrm{vv}}){\tilde{S}}_{\mathrm{hv}}^{*}\&rang;)}{-\&lang;{|{\tilde{S}}_{\mathrm{hh}}-{\tilde{S}}_{\mathrm{vv}}|}^2\&rang;+4\&lang;{\left|{\tilde{S}}_{\mathrm{hv}}\right|}^2\&rang;}\right)+\mathrm{\&pi;}]---\left(12\right)$

The computer capacity of arc tangent is [π, π] in the expression formula (12),

Obtain from the polarimetric synthetic aperture radar data;

The concrete steps that singular point is handled are as follows:

A) an angle change threshold T is set;

B) singularity of detection pixel, the order of pressing Row Column scans each pixel successively in entire image, 0 expression current pixel, 1,2,3,4 expressions are adjacent four pixels, if when the angle absolute difference of arbitrary pixel is greater than given threshold value T in the current pixel 0 and 1,4, judge that then current pixel point is a singular point;

C) if current pixel point is a non-singular point, then return step (b); Otherwise, enter step (d);

D) intermediate value with all pixel values in the given window size replaces current unusual pixel;

E) repeat (b)-(d) step 3 time, obtain the polarization orientation angle behind the singular point at last.

3, method according to claim 1 is characterized in that, adopts obtaining from the shade to the shape technology apart from the expression formula to gradient α in the radar angle measurement method:

$\mathrm{\&alpha;}=\mathrm{\&theta;}-\arcsin ((\sqrt{R^{\&prime;2}+4}-R^{\&prime;})/2)---\left(20\right)$

Wherein

$R^{\&prime;}=\frac{R(\mathrm{\&alpha;},\mathrm{\&beta;}){\cos }^2\mathrm{\&theta;}}{\sin \mathrm{\&theta;}\cos \mathrm{\&beta;}}---\left(21\right)$

The angle of depression θ of radar incident electromagnetic wave is obtained by the SAR data of input;

Again because polarization orientation angle ψ and terrain slope and radar visual angle have following relation

$\tan \mathrm{\&psi;}=\frac{\tan \mathrm{\&beta;}}{-\tan \mathrm{\&alpha;}\cos \mathrm{\&phi;}+\sin \mathrm{\&phi;}}---\left(13\right)$

Wherein, φ is the radar visual angle, and tan β is that the orientation is to gradient S _A, tan α is that distance is to gradient S _R,

From expression formula (13) and (20), (α β), thereby obtains the orientation to gradient S just can directly to calculate the surface slope angle _AWith the distance to gradient S _R

4, method according to claim 1 is characterized in that, definition weight quality figure from co-polarization signal response figure

$w_{i,j}=\frac{{\mathrm{\&sigma;}}_{\max }(0,{\mathrm{\&psi;}}^{\&prime;})-{\mathrm{\&sigma;}}_{\min }(0,{\mathrm{\&psi;}}^{\&prime;})}{{\mathrm{\&sigma;}}_{\max }(0,{\mathrm{\&psi;}}^{\&prime;})}---\left(27\right)$

σ wherein _Max(0, ψ ') and σ _Min(0, ψ ') remarked pixel (i, the j) maximum/minimum response of the co-polarization signal under linear polarization base (x=0) situation,

The calculation expression formula of co-polarization signal response figure is

$\mathrm{\&sigma;}(\mathrm{\&chi;},{\mathrm{\&psi;}}^{\&prime;})={\left[\begin{array}{c}1\\ \cos 2\mathrm{\&chi;}\cos 2{\mathrm{\&psi;}}^{\&prime;}\\ \cos 2\mathrm{\&chi;}\sin 2{\mathrm{\&psi;}}^{\&prime;}\\ \sin 2\mathrm{\&chi;}\end{array}\right]}^{T}K\left[\begin{array}{c}1\\ \cos 2\mathrm{\&chi;}\cos 2{\mathrm{\&psi;}}^{\&prime;}\\ \cos 2\mathrm{\&chi;}\sin 2{\mathrm{\&psi;}}^{\&prime;}\\ \sin 2\mathrm{\&chi;}\end{array}\right]---\left(28\right)$

Wherein, K is the Stokes matrix, and x is the ellipticity angle, and ψ ' is a polarization orientation angle.

5, method according to claim 1 is characterized in that adopting the least-squares algorithm of weighting to find the solution the discrete Poisson equation of relevant ground elevation.

Images