CN102183761A - Digital elevation model reconstruction method for space-borne interference synthetic aperture radar - Google Patents

Abstract

The invention relates to a digital elevation model reconstruction method for space-borne interference synthetic aperture radar. The method comprises the steps of: firstly, roughly sampling the original absolute interference phase diagram; and then, fitting to obtain the polynomials of the digital elevation model reconstruction three-dimensional coordinate at each pixel part of the rough sampling absolute interference phase diagram and the absolute interference phase; finally, obtaining the polynomial at each pixel part of the original absolute interference phase diagram by carrying out bilinear interpolation to the polynomial at each pixel part of the rough sampling absolute interference phase diagram; and substituting the absolute interference phase value at the pixel into the corresponding polynomial to obtain the corresponding three-dimensional coordinate. With less calculation, the invention can be finished in shorter time under the current computational ability; at the same time, the computational accuracy of the method can meet the demand of reconstructing the interference synthetic aperture radar digital elevation model.

Description

Spaceborne interference synthetic aperture radar digital elevation model method for reconstructing

Technical field

The invention belongs to the interleaving techniques field of microwave remote sensing and signal Processing, particularly a kind of interference synthetic aperture radar digital elevation model method for reconstructing.

Background technology

Interference synthetic aperture radar can obtain the digital elevation model with higher spatial resolution and vertical accuracy round-the-clock, round-the-clock.Utilize digital elevation model, can monitor the situation of change on earth top and ice and snow surface, carry out hazard predictions such as earthquake, volcanic explosion, landslide and great flood, for productions such as agricultural, forestry, fishery are offered help, for military activity provides the information support.In a word, the digital elevation model that utilizes interference synthetic aperture radar to obtain has important use value and application fields.

Early stage interference synthetic aperture radar is mainly used in mobile system, and bigger limitation is arranged on range of application and application space.The nineties in 20th century, interference synthetic aperture radar was successfully applied to spaceborne system along with development of technology, for the whole world provides abundant spaceborne interference of data of synthetic aperture radar.U.S. National Aeronautics and Space Administration and national image and Bureau of Surveying and Mapping united and carried out space shuttle radar mapping task on February 11st, 2000, test based on the interference synthetic aperture radar technology has specifically realized the three-dimensional mapping of global landform truly first, and the partial data of test and result have announced for each scientific research institution and handled and assessment.2007, Germany launched resolution and has reached 1 meter X-band Synthetic Aperture Radar satellite TerraSAR-X, and this is the highest civilian Synthetic Aperture Radar satellite of resolution up to now.Germany has launched second satellite about the same again on June 21st, 2010, two Synthetic Aperture Radar satellite will be formed the TanDEM-X system that double star forms into columns and carry out global topographic mapping, obtain global high accuracy number elevation model data, the expectation index of announcement will be higher than the space shuttle radar mapping task of the U.S..Germany has announced the repeat track interference of data of synthetic aperture radar of TerraSAR-X partly for global research institution in the internet and has used at present.

The digital elevation model reconstruction is the committed step in the interference of data of synthetic aperture radar processing procedure.The original echo data that the interference synthetic aperture radar system obtains through imaging, registration, interference, filtering, separate and obtain absolute interferometric phase after processing such as twining, digital elevation model is rebuild and is promptly utilized the orbit parameter of this absolute interferometric phase and radar system to come the pairing earth surface three-dimensional coordinate of each resolution element of inverting radar image.

The number with the radar image pixel consuming time that digital elevation model is rebuild is directly proportional.A vital task of spaceborne interference synthetic aperture radar system is to carry out whole world mapping, obtains global high accuracy number elevation model product.Germany space plan utilization TanDEM-X system of office finished global digital elevation model and measures in 3 years, estimate that the data volume of obtaining will be greater than 1300TB (the about 9.8TB of global metadata amount of u.s. space shuttle radar mapping task in 2000).In the face of magnanimity measured data like this, under the situation that guarantees precision, how improving data processing speed is a major issue that needs to be resolved hurrily.Therefore research is significant based on high precision, the fast digital elevation model method for reconstructing of spaceborne interference synthetic aperture radar system.

In the research of interference synthetic aperture radar digital elevation model method for reconstructing, the digital elevation model method for reconstructing that Goblirsch proposed under a kind of positive side-looking, the level land hypothesis in 1997, the three-dimensional coordinate error that this method reconstruction obtains is bigger; Giovanni Nico in 2002 based on distance ball, Doppler awl and the Phase Double curved surface analytic solution that digital elevation model is rebuild of having derived, the emission oblique distance that the derivation of these analytic solution is based on two width of cloth interference synthetic aperture radar images overlaps or the hypothesis of repeat track interferometry pattern; Marcus Studied the method for how obtaining elevation based on the product coordinate grid fast in 1998, this method can not obtain the three-dimensional coordinate of each resolution element of radar image; In 2005, Eineder carried out phase unwrapping and localization process simultaneously based on maximum Likelihood, and the interferometric phase image that merges under how much situations of different observations obtains digital elevation model.Above-mentioned digital elevation model method for reconstructing often is difficult to overcome processing speed and the contradiction of rebuilding precision when handling the magnanimity measured data.

Summary of the invention

The technical problem to be solved in the present invention provides interference synthetic aperture radar digital elevation model method for reconstructing, and this method calculated amount is little, can finish in the short period of time under existing computing power; The computational accuracy of this method can satisfy the requirement that interference synthetic aperture radar system digits elevation model is rebuild simultaneously.

The thinking of technical solution of the present invention is: at first, original absolute interferometric phase image is slightly sampled.Then, simulate the polynomial expression of each pixel place digital elevation model reconstruction of three-dimensional coordinate of the absolute interferometric phase image of thick sampling and absolute interferometric phase.At last, carry out bilinear interpolation by polynomial expression to each pixel place of the absolute interferometric phase image of thick sampling, obtain the polynomial expression at each pixel place of original absolute interferometric phase image, its corresponding polynomial expression of absolute interferometric phase value substitution at this pixel place can be obtained corresponding three-dimensional coordinate.

Technical solution of the present invention is:

Utilize the known original absolute interferometric phase image and the orbit parameter of radar system, finish following steps:

The first step: the absolute interferometric phase deviation map of deriving

To original absolute interferometric phase image according to the orientation to an interval k pixel, apart from (the k that slightly samples to an interval l pixel, the digital elevation model accuracy requirement decision that the value of l is obtained by final hope, the high more value of precision is more little, desirable [3,20] interval a certain round values), the absolute interferometric phase image of slightly being sampled, n phase deviation is superimposed to the absolute interferometric phase image of thick sampling respectively, derive the absolute interferometric phase deviation map of the n width of cloth, wherein the value of n is also determined by the digital elevation model accuracy requirement that final hope obtains, the high more value of precision is big more, a certain round values in desirable [4,10] are interval.When n is odd number, phase deviation V φ _i, i=1,2, the value of Ln is as follows:

$\left\{\begin{array}{cc}{\mathrm{V\φ}}_i=-i\·\frac{C\·2\mathrm{\π}\sqrt{{(f_a\·k)}^2+{(f_r\·l)}^2}}{n-1}& i=1,...\frac{n-1}{2},\\ V{\mathrm{\φ}}_i=0& i=\frac{n+1}{2},\\ V{\mathrm{\φ}}_i=-V{\mathrm{\φ}}_{i-\frac{n+1}{2}}& i=\frac{n+1}{2}+1,...n.\end{array}\right.$

When n is even number, phase deviation V φ _i, i=1,2, the value of Ln is as follows:

$\left\{\begin{array}{cc}{\mathrm{V\φ}}_i=-(2i-1)\·\frac{C\·2\mathrm{\π}\sqrt{{(f_a\·k)}^2+{(f_r\·l)}^2}}{2(n-1)}& i=1,...\frac{n}{2},\\ V{\mathrm{\φ}}_i=-V{\mathrm{\φ}}_{i-\frac{n}{2}}& i=\frac{n}{2}+1,...n.\end{array}\right.$

Wherein, C=1.5, f _a, f _rThe orientation that is respectively original absolute interferometric phase image to frequency and distance to frequency.

Second step: the digital elevation model of absolute interferometric phase deviation map is rebuild

Based on the absolute interferometric phase deviation map of each width of cloth, carry out digital elevation model in conjunction with oblique distance equation, Doppler equation and the interferometric phase equation of interference synthetic aperture radar system and rebuild, obtain n width of cloth digital elevation model reconstructed results.Digital elevation model is rebuild concrete grammar referring to Master's thesis " emulation and the applied research of the desirable amount of interference of space-based InSAR " (National University of Defense technology,, Wang Qingsong, the 21st page in 2009) 2.3.2 joint.

The 3rd step: fitting of a polynomial

Based on absolute interferometric phase deviation map of the aforementioned n width of cloth and corresponding n width of cloth digital elevation model reconstructed results, at each location of pixels (a of the absolute interferometric phase image of thick sampling _s, r _s) locate to simulate three the m rank polynomial expression that is independent variable with absolute interferometric phase value φ M ∈ [1, n-1] wherein, computational accuracy is as required chosen, polynomial expression

Be illustrated in location of pixels (a _s, r _s) x axial coordinate value in the digital elevation model locating to rebuild and the relation of absolute interferometric phase value φ,

Be illustrated in location of pixels (a _s, r _s) y axial coordinate value in the digital elevation model locating to rebuild and the relation of absolute interferometric phase value φ,

Be illustrated in location of pixels (a _s, r _s) z axial coordinate value in the digital elevation model locating to rebuild and the relation of absolute interferometric phase value φ.

The 4th step: polynomial interpolation

To polynomial expression

Carrying out bilinear interpolation can be in the hope of original each location of pixels of absolute interferometric phase image (a, three polynomial expression p ' that r) locate _{X (a, r)}(φ), p ' _{Y (a, r)}(φ), p ' _{Z (a, r)}(φ):

$p_{x(a,r)}^{\′}\left(\mathrm{\φ}\right)=\mathrm{bilinear}\left(\begin{array}{c}\left[\begin{array}{cc}{p}_{x(\mathrm{floor}\left(\frac{a}{k}\right),\mathrm{floor}\left(\frac{r}{l}\right))}\left(\mathrm{\φ}\right)& p_{x(\mathrm{floor}\left(\frac{a}{k}\right),\mathrm{floor}\left(\frac{r}{l}\right)+1)}\left(\mathrm{\φ}\right)\\ p_{x(\mathrm{floor}\left(\frac{a}{k}\right)+1,\mathrm{floor}\left(\frac{r}{l}\right))}\left(\mathrm{\φ}\right)& p_{x(\mathrm{floor}\left(\frac{a}{k}\right)+1,\mathrm{floor}\left(\frac{r}{l}\right)+1)}\left(\mathrm{\φ}\right)\end{array}\right]\end{array}\right)$

$p_{y(a,r)}^{\′}\left(\mathrm{\φ}\right)=\mathrm{bilinear}\left(\begin{array}{c}\left[\begin{array}{cc}{p}_{y(\mathrm{floor}\left(\frac{a}{k}\right),\mathrm{floor}\left(\frac{r}{l}\right))}\left(\mathrm{\φ}\right)& p_{y(\mathrm{floor}\left(\frac{a}{k}\right),\mathrm{floor}\left(\frac{r}{l}\right)+1)}\left(\mathrm{\φ}\right)\\ p_{y(\mathrm{floor}\left(\frac{a}{k}\right)+1,\mathrm{floor}\left(\frac{r}{l}\right))}\left(\mathrm{\φ}\right)& p_{y(\mathrm{floor}\left(\frac{a}{k}\right)+1,\mathrm{floor}\left(\frac{r}{l}\right)+1)}\left(\mathrm{\φ}\right)\end{array}\right]\end{array}\right)$

$p_{z(a,r)}^{\′}\left(\mathrm{\φ}\right)=\mathrm{bilinear}\left(\begin{array}{c}\left[\begin{array}{cc}{p}_{z(\mathrm{floor}\left(\frac{a}{k}\right),\mathrm{floor}\left(\frac{r}{l}\right))}\left(\mathrm{\φ}\right)& p_{z(\mathrm{floor}\left(\frac{a}{k}\right),\mathrm{floor}\left(\frac{r}{l}\right)+1)}\left(\mathrm{\φ}\right)\\ p_{z(\mathrm{floor}\left(\frac{a}{k}\right)+1,\mathrm{floor}\left(\frac{r}{l}\right))}\left(\mathrm{\φ}\right)& p_{z(\mathrm{floor}\left(\frac{a}{k}\right)+1,\mathrm{floor}\left(\frac{r}{l}\right)+1)}\left(\mathrm{\φ}\right)\end{array}\right]\end{array}\right)$

Wherein, the operation of bilinear () expression bilinear interpolation, floor () expression rounds operation to for a short time.

The 5th step: digital elevation model calculates

With original each location of pixels of absolute interferometric phase image (a, original absolute interferometric phase value φ (a, r) the substitution polynomial expression p ' that r) locate _{X (a, r)}(φ), p ' _{Y (a, r)}(φ), p ' _{Z (a, r)}(φ), can obtain the pairing terrain object point three-dimensional coordinate of this location of pixels (x, y, z), the digital elevation model reconstructed results that promptly finally obtains.

$\left\{\begin{array}{c}x=p_{x(a,r)}^{\′}\left(\mathrm{\φ}(a,r)\right)\\ y=p_{y(a,r)}^{\′}\left(\mathrm{\φ}(a,r)\right)\\ z=p_{z(a,r)}^{\′}\left(\mathrm{\φ}(a,r)\right)\end{array}\right.$

Adopt the present invention desirable following technique effect:

In the digital elevation model method for reconstructing of the present invention, the digital elevation model process of reconstruction of the absolute interferometric phase deviation map of several that derive is independently, can finish respectively, therefore digital elevation model method for reconstructing of the present invention extremely helps carrying out parallel computation, can give full play to the advantage of cluster computer, the design concurrent program significantly improves program run efficient.

Digital elevation model method for reconstructing of the present invention utilizes fitting of a polynomial greatly to reduce the number of times that iteration is found the solution the three dimensional non-linear system of equations, effectively reduces calculated amount.Therefore, adopt the present invention can finish digital elevation model efficiently and rebuild, be implemented under the existing computing power, finish digital elevation model fast and rebuild, precision can satisfy the requirement of interference synthetic aperture radar system mapping task simultaneously.

Accompanying drawing (table) explanation

Fig. 1 interference synthetic aperture radar digital elevation model provided by the invention method for reconstructing schematic flow sheet;

The original absolute interferometric phase image of Fig. 2 for using in the real data processing;

The digital elevation model figure that Fig. 3 obtains for traditional pointwise method for reconstructing;

Fig. 4 is the digital elevation model figure that utilizes the present invention to obtain.

Embodiment

Fig. 1 is an interference synthetic aperture radar digital elevation model method for reconstructing schematic flow sheet provided by the invention.Whole flow process was divided into for five steps.The first step, the absolute interferometric phase deviation map of deriving; In second step, the absolute interferometric phase deviation map of deriving is carried out digital elevation model rebuild; In the 3rd step, fitting of a polynomial simulates the polynomial expression of each pixel place digital elevation model reconstruction of three-dimensional coordinate of the absolute interferometric phase image of thick sampling and absolute interferometric phase value; The 4th step, polynomial interpolation, carrying out bilinear interpolation by the polynomial expression that match is obtained can be in the hope of the polynomial expression at each pixel place of original absolute interferometric phase image; In the 5th step, digital elevation model calculates, and can obtain final digital elevation model reconstructed results in the polynomial expression that the original absolute interferometric phase value substitution interpolation at each pixel place of original absolute interferometric phase image is obtained.

Fig. 2 is the original absolute interferometric phase image of input in the data processing experiment.Counting of these data is 4381 * 3386, about 9000 meters * 7000 meters of the scene domain that this width of cloth figure describes.

The digital elevation model figure that Fig. 3 obtains for traditional pointwise method for reconstructing, the processing time is 2149.016 seconds.

Fig. 4 is for utilizing the specific embodiment of the present invention at k=8, l=8, and n=4, the digital elevation model figure that obtains under the m=3 situation, the processing time is 231.012 seconds.Reconstructed results with respect to Fig. 3, Fig. 4 is 0.03745 meter at the error of coordinate root-mean-square value of x axle, error of coordinate root-mean-square value at the y axle is 0.18804 meter, error of coordinate root-mean-square value at the z axle is 0.17574 meter, with respect to original scene, rebuild the requirement that precision can satisfy interference synthetic aperture radar system mapping task.

Claims (1)

1. a spaceborne interference synthetic aperture radar digital elevation model method for reconstructing is characterized in that, utilizes the known original absolute interferometric phase image and the orbit parameter of radar system, finishes following steps:

The first step: the absolute interferometric phase deviation map of deriving;

Original absolute interferometric phase image is slightly sampled to an interval l pixel to an interval k pixel, distance according to the orientation, the digital elevation model accuracy requirement decision that the value of k, l is obtained by final hope, the high more value of precision is more little, the absolute interferometric phase image of slightly being sampled; N phase deviation is superimposed to the absolute interferometric phase image of thick sampling respectively, derive the absolute interferometric phase deviation map of the n width of cloth, wherein the value of n also determines that by the digital elevation model accuracy requirement that final hope obtains the high more value of precision is big more, a certain round values in desirable [4,10] are interval.When n is odd number, phase deviation V φ _i, i=1,2, the value of Ln is as follows:

When n is even number, phase deviation V φ _i, i=1,2, the value of Ln is as follows:

Wherein, C=1.5, f _a, f _rThe orientation that is respectively original absolute interferometric phase image to frequency and distance to frequency.

Second step: the digital elevation model of absolute interferometric phase deviation map is rebuild

Based on the absolute interferometric phase deviation map of each width of cloth, carry out digital elevation model in conjunction with oblique distance equation, Doppler equation and the interferometric phase equation of interference synthetic aperture radar system and rebuild, obtain n width of cloth digital elevation model reconstructed results.

The 3rd step: fitting of a polynomial

Based on absolute interferometric phase deviation map of the aforementioned n width of cloth and corresponding n width of cloth digital elevation model reconstructed results, at each location of pixels (a of the absolute interferometric phase image of thick sampling _s, r _s) locate to simulate three the m rank polynomial expression that is independent variable with absolute interferometric phase value φ

M ∈ [1, n-1] wherein, computational accuracy is as required chosen, polynomial expression

Be illustrated in location of pixels (a _s, r _s) x axial coordinate value in the digital elevation model locating to rebuild and the relation of absolute interferometric phase value φ,

Be illustrated in location of pixels (a _s, r _s) y axial coordinate value in the digital elevation model locating to rebuild and the relation of absolute interferometric phase value φ,

Be illustrated in location of pixels (a _s, r _s) z axial coordinate value in the digital elevation model locating to rebuild and the relation of absolute interferometric phase value φ.

The 4th step: polynomial interpolation

To polynomial expression

Carrying out bilinear interpolation can be in the hope of original each location of pixels of absolute interferometric phase image (a, three polynomial expression p ' that r) locate _{X (a, r)}(φ), p ' _{Y (a, r)}(φ), p ' _{Z (a, r)}(φ):

Wherein, the operation of bilinear () expression bilinear interpolation, floor () expression rounds operation to for a short time.

The 5th step: digital elevation model calculates

With original each location of pixels of absolute interferometric phase image (a, original absolute interferometric phase value φ (a, r) the substitution polynomial expression p ' that r) locate _{X (a, r)}(φ), p ' _{Y (a, r)}(φ), p ' _{Z (a, r)}(φ), can obtain the pairing terrain object point three-dimensional coordinate of this location of pixels (x, y, z), the digital elevation model reconstructed results that promptly finally obtains.

The digital elevation model accuracy requirement decision that the value of k, l is obtained by final hope, the high more value of precision is more little, desirable [3,20] interval a certain round values

The digital elevation model accuracy requirement decision that the value of n is also obtained by final hope, the high more value of precision is big more, a certain round values in desirable [4,10] are interval.

Images