CN103728604A - Broadband synthetic aperture radar sub-band interferometric data processing method - Google Patents

Abstract

Provided is a broadband synthetic aperture radar sub-band interferometric data processing method. The broadband synthetic aperture radar sub-band interferometric data processing method comprises the steps of (1) carrying out self-adaptive precise registration on main images and auxiliary images and resampling the auxiliary images, (2) optimizing sub-band frequency spectrum positions and broadband, (3) carrying out band-pass filtering on the main images and the auxiliary images, (4) generating sub-band interferograms and sub-band difference interferograms, (5) carrying out phase compensation on a reference plane and a DEM, and (6) carrying out stack processing on the sub-band difference interferograms. According to the broadband synthetic aperture radar sub-band interferometric data processing method, absolute radar interferometric phases can be obtained without phase unwrapping under most conditions, digital terrain maps or earth surface deformation maps can be directly obtained through transformation, and error diffusion caused by phase unwrapping mistakes can be eliminated. The broadband synthetic aperture radar sub-band interferometric data processing method is especially suitable for digital topographic mapping in high mapping regions, monitoring on landslides with large deformation velocity and monitoring on mine collapse.

Description

A kind of wideband synthetic aperture radar subband interference data disposal route

Technical field

The present invention relates to a kind of wideband synthetic aperture radar subband interference data disposal route, belong to interfering synthetic aperture radar technical field.It can obtain absolute radar interference phase place in most of the cases without phase unwrapping, thereby can directly obtain digital topography map or Ground Deformation figure by changing.

Background technology

Interfering synthetic aperture radar is the technology that a phase data according to high-resolution radar image is extracted terrain object three-dimensional spatial information.No matter utilize radar interference to carry out the measurement of landform digital elevation or Ground Deformation measurement, have the process of a phase unwrapping.Because the phase place of radar interference figure is all wound around by the angle of 2 π, so-called phase unwrapping, is exactly by the phase gradient between accumulation radar interference image primitive, thereby obtains the process of pixel absolute phase.This process is wanted to access correct result, and for topographical surveying, a necessary condition is that, in spatial domain, the absolute value of the phase differential between adjacent two pixels must not be greater than π; For the distortion measurement that utilizes the Time Series Analysis Method of storehouse differential interferometry figure, except that necessary condition must be met above, in time domain, the absolute value of the phase differential on adjacent two time points of same pixel also must not be greater than π.In Practical Project, especially when the topographical surveying of high Mountain area, or to having while carrying out interferometry compared with the Deformation Field in the landslide of large deformation velocity gradient and mine, the condition that the absolute value of phase differential above-mentioned must not be greater than π usually can not be met, thereby cause measurement result to occur very large error, sometimes even there will be without the situation of separating.According to radar interference theory, solve the problem of this phase unwrapping, to try every possible means exactly to improve interferometric phase and measure blur level.For topographical surveying, its blur level is defined as:

$H_{2\mathrm{\π}}=-\frac{\mathrm{\λr}}{{2B}_{\⊥}}\sin \left(\mathrm{\θ}\right)---\left(1\right)$

In above formula, λ is radar carrier wavelength, and r is the distance of terrain object and radar, and θ is radar incident angle, B _⊥for radar interference Space Baseline, H _{2 π}be elevation blur level, its represents interference phase difference corresponding ground elevation while being 2 π.

For distortion measurement, its blur level is defined as:

${\mathrm{\Δr}}_{2\mathrm{\π}}=\frac{\mathrm{\λ}}{2}---\left(2\right)$

In above formula, λ is radar carrier wavelength, △ r _{2 π}be deformation blur level, its represents differential interferometry phase differential corresponding Ground Deformation amount while being 2 π.

Obviously, from (1) and (2) formula, can find out, increase radar carrier wavelength and can increase elevation and the deformation blur level of interferometry, but this can increase hardware manufacturing difficulty and cost, and effect is not remarkable yet, because several times or ten increases doubly still can not meet actual demands of engineering.Present patent application is utilized the characteristic of wideband radar signal, in frequency domain, first radar image is split into several sub-band images, then corresponding radar sub-band images is interfered to computing, two subband interferograms of height are interfered again, at this moment simulating wavelength is an extended very large multiple, this multiple depends on the ratio of 2 times of original centre carrier frequency and subband center frequency, thereby can expand the wavelength of original centimetre-sized to meter level, the degree of difficulty of a difficult problem-phase unwrapping in radar interference technology will be eased, even can avoid this process of phase unwrapping.

Summary of the invention

1. object: the object of this invention is to provide a kind of wideband synthetic aperture radar subband interference data disposal route, it has overcome the deficiencies in the prior art, in most of the cases without phase unwrapping, can obtain absolute radar interference phase place, thereby can directly obtain digital topography map or Ground Deformation figure by changing, eliminate the error diffusion causing due to phase unwrapping mistake.

2. technical scheme: the present invention is a kind of wideband synthetic aperture radar subband interference data disposal route, and the method concrete steps are as follows:

Step 1: the self-adapting precision registration of major-minor image and the resampling of sub-picture

Because subband interference technique is to be mainly applied in the measurement of higher degree of high Mountain area or to have compared with the measurement of subsiding of the landslide of large deformation speed and mine, when major-minor image carries out registration, not only to consider radar interference system (imaging geometry, interfere how much, the skew of the position of the same target in ground radar orbit etc.) causing on major-minor image, also will consider the position skew on major-minor image that terrain object causes because of larger elevation or larger deformational displacement.In the present invention, sub-picture utilizes Sinc function to resample, and the displacement of Sinc function is determined by following formula:

Offset(i,j)=(Δi,Δj)=((Δi，Δj) _i,j) _System+((Δi,Δj) _i,j) _Object (3)

In above formula, offset (i, j) has represented to be positioned in sub-picture the pixel of the capable j row of i with respect to the side-play amount of same pixel position in master image, its side-play amount that to be system deviation and target itself cause because of elevation or deformation with.The displacement when neighborhood that offset (i, j) is positioned at the pixel of (i, j) in sub-picture s is got maximal value with the covariance function of the neighborhood of corresponding master image p pixel determines, sees (4) and (5) formula below

$c_{\mathrm{ps}}{(n,m)}_{i,j}=\underset{k=-N/2}{\overset{k=N/2}{\mathrm{\Σ}}}\underset{l=-M/2}{\overset{l=M/2}{\mathrm{\Σ}}}(p(k,l)-\tilde{p})\·(s(k+n,l+m)-\tilde{s})---\left(4\right)$

$\left\{\begin{array}{c}\mathrm{\Δi}=n_{\max }\\ \mathrm{\Δj}=m_{\max }\end{array}\right.---\left(5\right)$

In upper two formulas, p and s represent respectively major-minor image, c is illustrated in the covariance function of two windows of the capable M row of the N that classifies center as with the capable j of i in major-minor image, n and m represent respectively row and column side-play amount, k and l represent to participate in window the ranks position of correlation computations pixel, and (Δ i, Δ j) represents the ranks number when c obtains maximal value, i.e. (nmax, mmax).

Step 2: the optimization of subband spectrum position and bandwidth

In ideal conditions, we wish that subband bandwidth is wide as much as possible, can make like this radar image resolution loss not too big.On the other hand, we wish that upper and lower sub-band images does not have identical frequency content, otherwise these identical frequency contents can appear at the form of noise in secondary interferogram.And according to signal processing theory, ideal bandpass filter is not attainable, that is to say, input signal can only be subject to limited decay in the stopband of bandpass filter, and therefore, we wish that the position of upper and lower subband is at a distance of more far better.Obviously, two perfect conditions can not be met simultaneously.The selection of a compromise, is that the bandwidth bw of upper and lower subband equals 1/3rd of radar original signal bandwidth B, and upper and lower subband is positioned at f _c± f ₀position on,

$\left\{\begin{array}{c}\mathrm{bw}=\frac{B}{3}\\ f_0=\frac{B}{3}\end{array}\right.---\left(6\right)$

Under such optimization is selected, the ratio of simulating wavelength and carrier wavelength is:

$\frac{{\mathrm{\λ}}^{\′}}{\mathrm{\λ}}=\frac{{3f}_c}{2B}---\left(7\right)$

Generally:

$\left\{\begin{array}{c}{\mathrm{<figure-callout\; id="0"\; label="f"\; filenames="DEST\_PATH\_HDA0000465061860000012.png,DEST\_PATH\_HDA0000465061860000013.png"\; state="\{\{state\}\}">f</figure-callout>}}_0=\frac{B-\mathrm{bw}}{2}\\ {\mathrm{\&lambda;}}^{\&prime;}=\frac{f_c}{B-\mathrm{bw}}\&CenterDot;\mathrm{\&lambda;}\end{array}\right.---\left(8\right)$

In three formulas, B is radar pulse bandwidth in the above, and fc represents carrier frequency, and λ represents carrier wavelength, and bw represents upper and lower subband bandwidth, f ₀represent that upper and lower subband center frequency is for the side-play amount of carrier frequency, the simulating wavelength of λ ' expression subband differential interferometry figure.

Step 3: the bandpass filtering of major-minor image

Major-minor radar SLC image, by upper and lower subband band-pass filter, generates the complex image of subband up and down separately.The normalization centre frequency of upper subband bandpass filter is 0.333, and normalization bandwidth is 0.333; Lower subband bandpass filter normalization centre frequency is 0.666 ,-0.333, and normalization bandwidth is 0.333.In order to reduce the leakage of spectrum energy, the frequency spectrum of bandpass filter needs windowing, and window function is selected Caesar's window (Kaiser), and its spectrum expression formula is as follows:

$W(f;\mathrm{\&beta;},F)=\frac{I_0\left(\mathrm{\&beta;}\sqrt{1-{(2f/F)}^2}\right)}{I_0\left(\mathrm{\&beta;}\right)}-\frac{F}{2}<f<\frac{F}{2}---\left(9\right)$

In above formula, F is pass band width, I ₀be zero Bessel function, β is a constant controlling the wide and stopband attenuation tradeoff of window passband, conventionally gets between 3～6.

Fig. 2, Fig. 3 and Fig. 4 have shown respectively original radar SLC image and by the frequency spectrum in oblique distance (Slant Range) direction of the later sub-band images up and down of band-pass filter.

Step 4: subband interferogram and the map generalization of subband differential interferometry

The figure of subband up and down of master image interferes with the figure of subband up and down of corresponding sub-picture, subband interferogram and lower subband interferogram in generation.Upper and lower subband interferogram is interfered again, generate subband differential interferometry figure.Note P and S are respectively major-minor image, the subband differential interferometry figure that PS is them,

$\mathrm{PS}=(P_{\mathrm{up}}\&CenterDot;S_{\mathrm{up}}^{*})\&CenterDot;{(P_{\mathrm{low}}\&CenterDot;S_{\mathrm{low}}^{*})}^{*}---\left(10\right)$

In above formula, P _up, P _low, S _up, and S _lowrepresent respectively the figure of subband up and down being split off by major-minor image, PS represents the subband differential interferometry figure consisting of this fourth officer sub-band images.

Step 5: the phase compensation of reference surface and landform altitude (DEM)

In the present invention, the WGS84 ellipsoid model of the earth is as the reference surface of interferometry, if carry out differential interferometry measurement, the WGS84 spheroid surface of the earth also needs the DEM in the stack region of surveying, and the simulation interferometric phase that reference surface and landform altitude (DEM) cause will be removed from subband differential interferometry figure.For this reason, first need to calculate the distance R of each pixel and radar in major-minor image, it meets Nonlinear System of Equations below:

$\left\{\begin{array}{c}R=\sqrt{(\stackrel{\&RightArrow;}{R_s}-\stackrel{\&RightArrow;}{R_t})\&CenterDot;(\stackrel{\&RightArrow;}{R_s}-\stackrel{\&RightArrow;}{R_t})}\\ \frac{2}{\mathrm{\&lambda;R}}(\stackrel{\&RightArrow;}{V_s}-\stackrel{\&RightArrow;}{V_t})\&CenterDot;(\stackrel{\&RightArrow;}{R_s}-\stackrel{\&RightArrow;}{R_t})=0\\ \frac{x_t^2}{{(a+h)}^2}+\frac{y_t^2}{{(a+h)}^2}+\frac{z_t^2}{b^2}=1\end{array}\right.---\left(11\right)$

In above formula, R _sfor radar site vector, R _tfor pixel position vector, i.e. (x _t, y _t, z _t), V _sand V _tbe respectively radar and pixel target speed, h is pixel elevation, if compensate for reference face phase place is only got h=0, a and b are the major and minor axis of earth ellipsoid body Model.Calculate each pixel respectively in major-minor image with the distance R of radar after, the poor Δ R of the distance by pixel in major-minor, just can calculate reference surface in subband differential interferometry figure and the analogue phase of landform altitude:

$\mathrm{\&phi;}=-\frac{4\mathrm{\&pi;}}{{\mathrm{\&lambda;}}^{\&prime;}}\&CenterDot;\mathrm{\&Delta;R}---\left(12\right)$

Wavelength in above formula should be the simulating wavelength after subband is interfered.Fig. 5 is the schematic diagram of the main process of the formation since the radar image after two registrations to subband differential interferometry figure.Utilize

(1) formula and

(2) expressed elevation and the deformation blur level of formula just can be converted into earth's surface elevation or Surface Deformation Field by subband differential interferometry figure.

Step 6: the storehouse of subband differential interferometry figure is processed

A) generation of digital elevation figure (DEM)

Although the elevation blur level of subband differential interferometry figure has been enhanced doubly a lot, also the variance of the error of interferometric phase has also been expanded and has been close to identical multiple simultaneously, the phase place of individual subband differential interferometry figure can contain larger noise.Because the phase place of interferogram is proportional to elevation by the factor of phase ambiguity, therefore multiple subband differential interferometry figure can be formed to a storehouse, then by least square method, estimate the elevation of each pixel:

$\underset{h}{\mathrm{Min}}\mathrm{\&Sigma;}{(h-\frac{\mathrm{\&Delta;}{\mathrm{\&phi;}}_i}{2\mathrm{\&pi;}}{H}_{2\mathrm{\&pi;},i})}^2---\left(13\right)$

In above formula, Δ φ i represents pixel i phase place in a certain interferogram in storehouse, H _{2 π, i}the elevation mould that represents this point

Paste degree, h represents the estimated value of this point height.

B) measurement of Surface Deformation Field

Because deformation (displacement) blur level of subband differential interferometry figure is very large, therefore while being converted to deformation quantity by differential phase hardly with carrying out phase unwrapping, but in order to eliminate atmospheric effect, vertical error and other various noises, still be necessary to utilize the time series analysis of one group of subband differential interferometry figure (storehouse), particularly short baseline set method, tries to achieve accurate Surface Deformation Field and dynamic change thereof.It is identical substantially that the Time Series Analysis Method of subband differential interferometry figure is interfered time Sequence Analysis Method with traditional full bandwidth data, unique different be the choosing method of high coherent point.In the present invention, utilized two criterions to determine candidate's high coherent point:

$\left\{\begin{array}{c}\frac{\mathrm{\&sigma;}}{\mathrm{\&mu;}}<c\\ \mathrm{SCR}>d\end{array}\right.---\left(14\right)$

In above formula, μ and σ are average and the variances of pixel intensity level in storehouse time range, c is the threshold value that the ratio of variance and average is set, and it can weigh this pixel whether stablize of amplitude in storehouse time range substantially, and generally value is 0.2 left and right.SCR is the signal to noise ratio (Signal to Clutter Ratio) of pixel, and in subband differential interferometry figure, the signal to noise ratio that can carry out the pixel of time series analysis must be greater than a threshold value d, general value 25dB.Chosen after candidate point, utilized general period map method or short baseline set method can obtain the speed of Surface Deformation Field, linear and non-linear deformation quantity.If by artificial reverberator, the sudden change of large scale (tens of centimetres even meter level) and the measurement of non-linear deformation can not be restricted in the present invention.

3. advantage and effect: the invention provides a kind of wideband synthetic aperture radar subband interference data disposal route, its advantage is:

(1) this method utilizes the subband signal of wideband radar to interfere, become hundred times of ground to increase analog carrier wavelength and interferometry elevation blur level, in high Mountain area with a varied topography, need not very easily cause that the wrong solution process of twining can obtain radar interference absolute phase, thereby obtain the digital elevation on earth's surface.

(2) this method is utilized the subband signal differential interferometry of wideband radar, can directly obtain the even deformation of meter level of tens of centimetres, earth's surface, has greatly expanded the measurement range of traditional differential interferometry radar.

(3) this method is utilized the sequence subband signal differential interferometry of wideband radar, the non-linear deformation that can obtain earth's surface large scale, and this is very effective in the deformation monitoring of landslide.

(4) due to this method, construct the subband signal of two different center frequency, thereby can as double-frequency GPS receiver, measure ionospheric electron density and change, eliminated the phase error that this variation causes, improved interferometry precision.

Accompanying drawing explanation

Fig. 1. the frequency band of wideband radar SLC image and division subband schematic diagram

In Fig. 1, f _cfor radar carrier frequency, f ₀for the side-play amount of subband center frequency with respect to radar carrier frequency, bw represents bandwidth.

Fig. 2. wideband radar SLC image distance is to frequency spectrum

Fig. 3. upper sub-band images distance is to frequency spectrum

Fig. 4. lower sub-band images distance is to frequency spectrum

Fig. 5. subband is interfered process flow diagram

Fig. 6. the TerraSAR-X radar image that tree level ground is come down regional

Fig. 7. No. 8 and three kinds of measuring method result comparison diagrams of No. 9 reverberator slide displacement amounts on the sliding mass of tree level ground

Fig. 8. FB(flow block) of the present invention

Embodiment

Landslide, tree level ground, dam area, the Yangtze River Gorges deformation monitoring of take is example, and the concrete operation step of the present invention in practical engineering application is described.See Fig. 8, the present invention is a kind of wideband synthetic aperture radar subband interference data disposal route, and the method concrete steps are as follows:

Step 1: the self-adapting precision registration of data decimation and major-minor image

1. subband interference technique applicability is considered

According to the topography and geomorphology of engineering purpose and workspace, ground vegetation level of coverage and Ground Deformation situation judge whether to enable subband interference technique.For take, to generate DEM be object, is only applicable to high Mountain area, for deformation monitoring, is only applicable to have the area of large deformation speed.No matter which kind of object, all requires ground monitored (or mapping) target can have higher radar coherence.For the deformation monitoring of the high coherent point by time Sequence Analysis Method, require ground monitored (or mapping) target can there is higher radar signal to noise ratio (SCR>25dB), select if desired and install artificial reverberator.Landslide, tree level ground, Three Gorges landform changes greatly, and crops and fruit tree cover throughout the year, and landslide unstable region slippage is very large, and conventional interference radar cannot be monitored out the deformation data on this landslide, therefore should adopt, artificial reverberator and subband differential interference method is installed.On the sliding mass of tree level ground, 14 artificial reverberators have been installed, have separately been had 4 to be arranged on outside sliding mass, as the benchmark of measuring and verifying.

2. data decimation

Choose the spaceborne wideband synthetic aperture radar SLC image that several cover workspace, more than bandwidth requirement 100MHz.For this purpose, the 29 width SLC data on the April 9,22 days to 2010 April in 2009 of 3 meters of resolution of the spaceborne X-band of German TerraSAR-X have been chosen.30 kilometers of data fabric widths, long 60 kilometers, bandwidth 150MHz, data cover time span 1 year, it is 11 days that satellite repeats the cycle of passing by.Radar general data parameter and scanning the date in Table 1 and table 2 shown in:

Table 1: select TerraSAR-X radar major parameter table

Carrier frequency	9.6500000+9Hz
		Carrier wavelength	3.1cm
Pulse bandwidth	1.5000000e+8Hz
		Incident angle	26.2542°
Distance is to Pixel size	0.909403m
		Orientation is to Pixel size	1.965121m
Distance is to sampling rate	1.6482919+8Hz
		Radar scanning pattern	Band
Data type	Haplopia plural number

Table 2: radar data date table

20090422	20090503	20090514	20090525	20090605	20090616	20090627
							20090708	20090719	20090730	20090810	20090821	20090901	20090912
20090923	20091015	20091026	20091117	20091128	20091209	20091220
							20091231	20100111	20100122	20100202	20100213	20100307	20100409

3. self-adapting precision registration

Choosing piece image is master image, and as with reference to image, it should be positioned at the center of time domain (time basis collection) and spatial domain (Space Baseline collection) in principle.Because the present invention goes for large scale distortion measurement, in tree level ground landslide monitoring again by means of artificial reverberator, therefore choosing of master image is not subject to the restriction chatted above.Convenient for time series analysis, choose first image, the data on January 7th, are 1 master image.Take (1) formula as basis, utilize (2) formula and (3) formula to calculate in all the other each sub-pictures all pixels with respect to the relativity shift vector of master image, with the Sinc function of displacement, each sub-picture is carried out to self-adaption two-dimensional resampling, to reach the object of precision registration, precision is 1/10th pixels.

Offset(i,j)=(Δi，Δj)=((Δi，Δj) _i,j) _System+((Δi,Δj) _i,j) _Object (1)

$c_{\mathrm{ps}}{(n,m)}_{i,j}=\underset{k=-N/2}{\overset{k=N/2}{\mathrm{\&Sigma;}}}\underset{l=-M/2}{\overset{l=M/2}{\mathrm{\&Sigma;}}}(p(k,l)-\tilde{p})\&CenterDot;(s(k+n,l+m)-\tilde{s})---\left(2\right)$

$\left\{\begin{array}{c}\mathrm{\&Delta;i}=n_{\max }\\ \mathrm{\&Delta;j}=m_{\max }\end{array}\right.---\left(3\right)$

In formula (1) first departure that represents to be arranged in master image the position of the pixel of the capable j row of i and the same place of sub-picture, this departure is a bivector, in formula (1) second.Bias vector is that the component that the own deformation of terrain object that in the component that produced by system and master image, this pixel is corresponding produces forms, and is respectively in formula (1) shown in third and fourth.

In the situation that the parameters such as satellite orbit number and radar imagery time are known, can first calculate the offset component that system produces.Take this systematic component as basis, can roughly determine pixel position corresponding in sub-picture.Like this, for each pixel (i, j) in master image, all can find the window of a suitable size centered by this pixel, such as 512 row and 512 row, can also determine the position of this window corresponding in sub-picture simultaneously.Calculate the covariance function of these two windows, as the formula (2).M, the ranks number that N is window, p and s represent respectively major-minor image.M while making covariance function obtain maximal value and n are the side-play amount of pixel (i, j) in sub-picture in master image, as the formula (3).

Step 2: Domain Design-subband spectrum position of bandpass filter and the optimization of bandwidth

At radar carrier frequency f _cin the known situation of radar pulse signal bandwidth B, determine the centre frequency f of bandpass filter ₀, the transport function of bandwidth bw and wave filter.First, it is far away that upper and lower two bandpass filter will apart be tried one's best, i.e. a least significant end in radar bandwidth spectrum, and another can increase to greatest extent like this bandpasstilter stopband decay, thereby reduce the noise in subband differential interferometry figure at most significant end.The simulating wavelength of large young pathbreaker's determinant band differential interferometry figure of centre frequency, also will determine the size of simulation blur level.Generally get bw=B/3.Because subband bandwidth is dwindled three times, obviously resolution can increase three times.After bw determines, f ₀=(B-bw)/2, the simulating wavelength of subband differential interferometry figure is k times of radar carrier wavelength, k=f _c/ (B-bw).According to the landform of workspace or deformation data maximal value, whether the elevation blur level or the deformation blur level that check subband interference be suitable, otherwise, can adjust by changing the size of bandpass filter centre frequency and bandwidth the size of k, thereby determine suitable blur level, to avoid the phase unwrapping of subband differential interferometry figure.In Domain Design, often by frequency spectrum to bandwidth B normalization, obtain like this optimal design value of bandpass filter below:

$\left\{\begin{array}{c}\mathrm{bw}=\frac{1}{3}\\ f_{\mathrm{up},0}=-f_{\mathrm{low},0}=\frac{1}{3}\\ k=\frac{{3f}_c}{2B}\end{array}\right.---\left(4\right)$

The window function of bandpass filter, is also transport function here, selects the Caesar's window shown in (5) formula, brings the parameter in formula (4) into formula (5), gets β=6.

$W(f;\mathrm{\&beta;},F)=\frac{I_0\left(\mathrm{\&beta;}\sqrt{1-{(2f/F)}^2}\right)}{I_0\left(\mathrm{\&beta;}\right)}-\frac{F}{2}<f<\frac{F}{2}---\left(5\right)$

In above formula, F is pass band width, I ₀be zero Bessel function, β controls the wide and stopband attenuation compromise of window passband

A constant of relation, gets between 3～6 conventionally

Step 3: the bandpass filtering of major-minor image and the generation of sub-band images

The generation of sub-band images is all carried out in frequency domain.First by image (p, s after all 29 registrations ₁..., s ₂₈) by the quick Fourier transforms (2D-FFT) of two dimension, convert frequency domain figure picture (P, S to ₁..., S ₂₈).Two band-pass filter up and down that these 29 frequency domain figure pictures design by step 4 respectively, generate upper sub-band images and lower sub-band images (P in frequency domain separately _up, P _low, S _{1_up}, S _{1_low}..., S _{28_up}, S _{28_low}).Before filtering, sub-band images will be also spectrum shift by demodulation, will be with logical spectrum shift on base band, and making the analog carrier frequency of all upper sub-band images is f _c+ f ₀, the analog carrier frequency of all lower sub-band images is f _c-f ₀, and the frequency spectrum of all sub-band images is all centered by its analog carrier frequency, to guarantee the correct of the same and result of frequency spectrum that follow-up signal processes.And then utilize the fast contrary Fourier transforms (2D-IFFT) of two dimension that these 58 frequency domain figures are looked like to convert to 58 spatial domain image (p _up, p _low, s _{1_up}, s _{1_low}..., s _{28_up}, s _{28_low}).

Fig. 1 is frequency band and the division subband schematic diagram of wideband radar SLC image, Fig. 2, and Fig. 3 and Fig. 4 have shown respectively original radar SLC image and have passed through the frequency spectrum in the later oblique distance of sub-band images up and down (Slant Range) direction of band-pass filter.

Step 4: subband interferogram and the map generalization of subband differential interferometry

On April 22nd, 2009 image of take is master image, and all the other are all the mode of sub-picture, generates 28 upper subband interferograms and 28 lower subband interferograms:

$p_{\mathrm{up}}\&CenterDot;s_{i\_\mathrm{up}}^{*},i=1,......,28---\left(6\right)$

$p_{\mathrm{low}}\&CenterDot;s_{i\_\mathrm{low}}^{*}.i=1,......,28---\left(7\right)$

$(p_{\mathrm{up}}\&CenterDot;s_{i\_\mathrm{up}}^{*})\&CenterDot;{(p_{\mathrm{low}}\&CenterDot;s_{i\_\mathrm{low}}^{*})}^{*},i=1,......,28---\left(8\right)$

As shown in (8) formula, each is interfered producing a subband differential interferometry figure, and it interferes interference image and two interference images that above sub-band images produces that two right lower sub-band images produce again to interfere and produce by this.Due to the bandwidth of subband differential interferometry figure now, only have 1/3rd original (in the situations of optimal design), pixel resolution has reduced by three times, therefore doing while again interfering, should do the processing of looking of 3x3 at least, to reduce speckle noise more.

Step 5: the phase compensation of reference surface and landform altitude

Phase compensation in the phase compensation procedure of reference surface and landform altitude and traditional differential interferometry substantially just as.If only do the generation of the digital elevation (DEM) of landform, subband differential interferometry figure only need make reference the phase compensation of face, if do Ground Deformation monitoring, except making reference the phase compensation of face, also needs to do the phase compensation of landform altitude.For the orbit data of the satellite with as radar carrier match, the surface of choosing earth WGS84 ellipsoid model is as with reference to face, the simulation interferometric phase that spheroid surface and the landform altitude above it produce is according to satellite orbit data, radar imagery phase data and (9) formula and (10) formula are calculated, different from traditional differential interferometry, at this moment guinea pig carrier wavelength is k times of primary carrier wavelength, and the value of k is shown in formula (11).From subband differential interferometry figure, deduct the simulation interferometric phase that this is produced by reference surface and landform altitude.If only estimate landform altitude, making the h in (9) formula is zero.

$\left\{\begin{array}{c}R=\sqrt{(\stackrel{\&RightArrow;}{R_s}-\stackrel{\&RightArrow;}{R_t})\&CenterDot;(\stackrel{\&RightArrow;}{R_s}-\stackrel{\&RightArrow;}{R_t})}\\ \frac{2}{\mathrm{\&lambda;R}}(\stackrel{\&RightArrow;}{V_s}-\stackrel{\&RightArrow;}{V_t})\&CenterDot;(\stackrel{\&RightArrow;}{R_s}-\stackrel{\&RightArrow;}{R_t})=0\\ \frac{x_t^2}{{(a+h)}^2}+\frac{y_t^2}{{(a+h)}^2}+\frac{z_t^2}{b^2}=1\end{array}\right.---\left(9\right)$

In above formula, R _sfor radar site vector, R _tfor pixel position vector, i.e. (x _t, y _t, z _t), V _sand V _tbe respectively radar and pixel target speed, h is pixel elevation, if compensate for reference face phase place is only got h=0, a and b are the major and minor axis of earth ellipsoid body Model.Calculate each pixel respectively in major-minor image with the distance R of radar after, the poor Δ R of the distance by pixel in major-minor, just can calculate reference surface in subband differential interferometry figure and the analogue phase of landform altitude:

$\mathrm{\&phi;}=-\frac{4\mathrm{\&pi;}}{{\mathrm{\&lambda;}}^{\&prime;}}\&CenterDot;\mathrm{\&Delta;R}---\left(10\right)$

Wavelength in above formula should be the simulating wavelength after subband is interfered:

$k=\frac{{3f}_c}{2B}---\left(11\right)$

Fig. 5 is the main process schematic diagram of the formation since the radar image after two registrations to subband differential interferometry figure.

Step 6: the storehouse of subband differential interferometry figure is processed

1. the filtering of subband differential interferometry figure

Due to the Gaussian characteristics of subband differential interferometry figure phase place, should adopt 2-d gaussian filters device to carry out filtering to the phase place of subband differential interferometry figure, eliminate Gaussian noise.

2. the measurement of the processing of the storehouse of subband differential interferometry figure and Surface Deformation Field

Because deformation (displacement) blur level of subband differential interferometry figure is very large, therefore while being converted to deformation quantity by differential phase hardly with carrying out phase unwrapping, but in order to eliminate atmospheric effect, vertical error and other various noises, still be necessary to utilize the time series analysis of one group of subband differential interferometry figure (storehouse), particularly short baseline set method, tries to achieve accurate Surface Deformation Field and dynamic change thereof.It is identical substantially that the Time Series Analysis Method of subband differential interferometry figure is interfered time Sequence Analysis Method with traditional full bandwidth data, unique different be the choosing method of high coherent point.In the present invention, utilized two criterions to determine candidate's high coherent point:

$\left\{\begin{array}{c}\frac{\mathrm{\&sigma;}}{\mathrm{\&mu;}}<c\\ \mathrm{SCR}>d\end{array}\right.---\left(12\right)$

In above formula, μ and σ are average and the variances of pixel intensity level in storehouse time range, c is the threshold value that the ratio of variance and average is set, and it can weigh this pixel whether stablize of amplitude in storehouse time range substantially, and generally value is 0.2 left and right.SCR is the signal to noise ratio (Signal to Clutter Ratio) of pixel, and in subband differential interferometry figure, the signal to noise ratio that can carry out the pixel of time series analysis must be greater than a threshold value d, general value 25dB.Chosen after candidate point, utilized general period map method or short baseline set method can obtain the speed of Surface Deformation Field, linear and non-linear deformation quantity.If by artificial reverberator, the sudden change of large scale (tens of centimetres even meter level) and the measurement of non-linear deformation can not be restricted in the present invention.

(12) formula of utilization, can determine the position of all 18 artificial reverberators in each subband differential interferometry figure, reads the phase place of these 18 artificial reverberators in every subband differential interferometry figure, and deducts the phase place of a reference point.In landslide, tree level ground, No. 6 artificial reverberator is chosen to be reference point.Obtain after the phase place of each artificial reverberator in each subband differential interferometry figure, by formula (13), can calculate each artificial reverberator at corresponding interferogram the displacement in the radar ray direction in the time period:

$\mathrm{\&rho;}=\frac{{\mathrm{\&lambda;}}^{\&prime;}}{4\mathrm{\&pi;}}\&CenterDot;\mathrm{\&Delta;\&phi;}---\left(13\right)$

Fig. 6 is the TerraSAR-X radar degree figure on landslide, tree level ground, and 14 artificial reverberators and the artificial reverberator outside 4 sliding masses on sliding mass are high-visible.In order to verify measurement accuracy of the present invention, chosen the subband differential interferometry result of No. 8 and No. 9 reverberators here and the synchro measure result of ground level and GPS is made comparisons.No. 8 reverberator is positioned at the outer stabilized zone of sliding mass, and No. 9 reverberator is positioned at sliding mass compared with large deformation district.Fig. 7 is the comparison diagram of result of three kinds of measuring methods of two reverberators.Comparative result demonstration, three kinds of measurement results are very identical, will particularly point out here, and during May to July, the displacement of No. 9 reverberators reaches nearly 40 centimetres, and this cannot measure by traditional radar interference method.

3. the subsequent treatment of subband differential interferometry figure

The solution of subband differential interferometry figure twines, the conversion of phase place-elevation, phase place-deformation conversion, the time series analysis of subband differential interferometry figure storehouse, the processing procedures such as geocoding are completely identical with traditional full bandwidth radar interference data processing, and the data processing of these several parts is not within the scope of the present invention.What make an exception is only to make the carrier wavelength in (14) and (15) formula blur level expression formula into subband to interfere simulating wavelength.For topographical surveying, its blur level is defined as:

$H_{2\mathrm{\&pi;}}=-\frac{\mathrm{\&lambda;r}}{{2B}_{\&perp;}}\sin \left(\mathrm{\&theta;}\right)---\left(14\right)$

In above formula, λ is radar carrier wavelength, and r is the distance of terrain object and radar, and θ is radar incident angle, B _⊥for radar interference Space Baseline, H _{2 π}be elevation blur level, its represents interference phase difference corresponding ground elevation while being 2 π.For distortion measurement, its blur level is defined as:

${\mathrm{\&Delta;r}}_{2\mathrm{\&pi;}}=\frac{\mathrm{\&lambda;}}{2}---\left(15\right)$

In above formula, λ is radar carrier wavelength, Δ r _{2 π}be deformation blur level, its represents differential interferometry phase differential corresponding Ground Deformation amount while being 2 π.

Claims (1)

1. a wideband synthetic aperture radar subband interference data disposal route, is characterized in that: the method concrete steps are as follows:

Step 1: the self-adapting precision registration of major-minor image and the resampling of sub-picture

Because subband interference technique is to be mainly applied in the measurement of higher degree of high Mountain area or to have compared with the measurement of subsiding of the landslide of large deformation speed and mine, when major-minor image carries out registration, not only to consider the skew of the same target in ground that radar interference system the causes position on major-minor image, also will consider the position skew on major-minor image that terrain object causes because of larger elevation or larger deformational displacement; Here, sub-picture utilizes Sinc function to resample, and the displacement of Sinc function is determined by following formula:

Offset(i,j)=(△i,△j)=((△i,△j) _i,j) _System+((△i,△j) _i,j) _Object (3)

In above formula, offset (i, j) has represented to be positioned in sub-picture the pixel of the capable j row of i with respect to the side-play amount of same pixel position in master image, its side-play amount that to be system deviation and target itself cause because of elevation or deformation with; The displacement when neighborhood that offset (i, j) is positioned at the pixel of (i, j) in sub-picture s is got maximal value with the covariance function of the neighborhood of corresponding master image p pixel determines, sees (4) and (5) formula below

$c_{\mathrm{ps}}{(n,m)}_{i,j}=\underset{k=-N/2}{\overset{k=N/2}{\mathrm{\&Sigma;}}}\underset{l=-M/2}{\overset{l=M/2}{\mathrm{\&Sigma;}}}(p(k,l)-\tilde{p})\&CenterDot;(s(k+n,l+m)-\tilde{s})---\left(4\right)$

$\left\{\begin{array}{c}\mathrm{\&Delta;i}=n_{\max }\\ \mathrm{\&Delta;j}=m_{\max }\end{array}\right.---\left(5\right)$

In upper two formulas, p and s represent respectively major-minor image, c is illustrated in the covariance function of two windows of the capable M row of the N that classifies center as with the capable j of i in major-minor image, n and m represent respectively row and column side-play amount, k and l represent to participate in window the ranks position of correlation computations pixel, and (△ i, △ j) represents the ranks number when c obtains maximal value, i.e. (nmax, mmax);

Step 2: the optimization of subband spectrum position and bandwidth

In ideal conditions, wish that subband bandwidth is wide as much as possible, can make like this radar image resolution loss not too big; On the other hand, wish that upper and lower sub-band images does not have identical frequency content, otherwise these identical frequency contents can appear at the form of noise in secondary interferogram; And according to signal processing theory, ideal bandpass filter is not attainable, that is to say, input signal can only be subject to limited decay in the stopband of bandpass filter, therefore, wishes that the position of upper and lower subband is at a distance of more far better; Obviously, two perfect conditions can not be met simultaneously, and the selection of a compromise is that the bandwidth bw of upper and lower subband equals 1/3rd of radar original signal bandwidth B, and upper and lower subband is positioned at f _c± f ₀position on,

$\left\{\begin{array}{c}\mathrm{bw}=\frac{B}{3}\\ f_0=\frac{B}{3}\end{array}\right.---\left(6\right)$

Under such optimization is selected, the ratio of simulating wavelength and carrier wavelength is:

$\frac{{\mathrm{\&lambda;}}^{\&prime;}}{\mathrm{\&lambda;}}=\frac{3f_c}{2B}---\left(7\right)$

Generally:

$\left\{\begin{array}{c}{f}_0=\frac{B-\mathrm{bw}}{2}\\ {\mathrm{\&lambda;}}^{\&prime;}=\frac{f_c}{B-\mathrm{bw}}\&CenterDot;\mathrm{\&lambda;}\end{array}\right.---\left(8\right)$

In three formulas, B is radar pulse bandwidth in the above, and fc represents carrier frequency, and λ represents carrier wavelength, and bw represents upper and lower subband bandwidth, f ₀represent that upper and lower subband center frequency is for the side-play amount of carrier frequency, the simulating wavelength of λ ' expression subband differential interferometry figure;

Step 3: the bandpass filtering of major-minor image

Major-minor radar SLC image, by upper and lower subband band-pass filter, generates the complex image of subband up and down separately; The normalization centre frequency of upper subband bandpass filter is 0.333, and normalization bandwidth is 0.333; Lower subband bandpass filter normalization centre frequency is 0.666 ,-0.333, and normalization bandwidth is 0.333; In order to reduce the leakage of spectrum energy, the frequency spectrum of bandpass filter needs windowing, and it is Kaiser that window function is selected Caesar's window, and its spectrum expression formula is as follows:

$W(f;\mathrm{\&beta;},F)=\frac{I_0\left(\mathrm{\&beta;}\sqrt{1-{(2f/F)}^2}\right)}{I_0\left(\mathrm{\&beta;}\right)}-\frac{F}{2}<f<\frac{F}{2}---\left(9\right)$

In above formula, F is pass band width, I ₀be zero Bessel function, β is a constant controlling the wide and stopband attenuation tradeoff of window passband, conventionally gets between 3～6;

Step 4: subband interferogram and the map generalization of subband differential interferometry

The figure of subband up and down of master image interferes with the figure of subband up and down of corresponding sub-picture, subband interferogram and lower subband interferogram in generation; Upper and lower subband interferogram is interfered again, generate subband differential interferometry figure; Note P and S are respectively major-minor image, the subband differential interferometry figure that PS is them,

$\mathrm{PS}=(P_{\mathrm{up}}\&CenterDot;S_{\mathrm{up}}^{*})\&CenterDot;{(P_{\mathrm{low}}\&CenterDot;S_{\mathrm{low}}^{*})}^{*}---\left(10\right)$

In above formula, P _up, P _low, S _up, and S _lowrepresent respectively the figure of subband up and down being split off by major-minor image, PS represents the subband differential interferometry figure consisting of this fourth officer sub-band images;

Step 5: the phase compensation of reference surface and landform altitude (DEM)

The WGS84 ellipsoid model of the earth is as the reference surface of interferometry, if carry out differential interferometry measurement, the WGS84 spheroid surface of the earth also needs the DEM in the stack region of surveying, the simulation interferometric phase that reference surface and landform altitude (DEM) cause will be removed from subband differential interferometry figure, for this reason, first need to calculate the distance R of each pixel and radar in major-minor image, it meets Nonlinear System of Equations below:

$\left\{\begin{array}{c}R=\sqrt{(\stackrel{\&RightArrow;}{R_s}-{\stackrel{\&RightArrow;}{R}}_t)\&CenterDot;(\stackrel{\&RightArrow;}{R_s}-{\stackrel{\&RightArrow;}{R}}_t)}\\ \frac{2}{\mathrm{\&lambda;R}}(\stackrel{\&RightArrow;}{V_s}-{\stackrel{\&RightArrow;}{V}}_t)\&CenterDot;(\stackrel{\&RightArrow;}{R_s}-{\stackrel{\&RightArrow;}{R}}_t)=0\\ \frac{x_t^2}{{(a+h)}^2}+\frac{y_t^2}{{(a+h)}^2}+\frac{z_t^2}{b^2}=1\end{array}\right.---\left(11\right)$

In above formula, R _sfor radar site vector, R _tfor pixel position vector, i.e. (x _t, y _t, z _t), V _sand V _tbe respectively radar and pixel target speed, h is pixel elevation, if compensate for reference face phase place is only got h=0, a and b are the major and minor axis of earth ellipsoid body Model; Calculate each pixel respectively in major-minor image with the distance R of radar after, the poor △ R of the distance by pixel in major-minor, just can calculate reference surface in subband differential interferometry figure and the analogue phase of landform altitude:

$\mathrm{\&phi;}=-\frac{4\mathrm{\&pi;}}{{\mathrm{\&lambda;}}^{\&prime;}}\&CenterDot;\mathrm{\&Delta;R}---\left(12\right)$

Wavelength in above formula should be the simulating wavelength after subband is interfered, and utilizes

(1) formula and

(2) expressed elevation and the deformation blur level of formula just can be converted into earth's surface elevation or Surface Deformation Field by subband differential interferometry figure;

Step 6: the storehouse of subband differential interferometry figure is processed

A) digital elevation figure is the generation of DEM

Although the elevation blur level of subband differential interferometry figure has been enhanced doubly a lot, also the variance of the error of interferometric phase has also been expanded and has been close to identical multiple simultaneously, the phase place of individual subband differential interferometry figure can contain larger noise; Because the phase place of interferogram is proportional to elevation by the factor of phase ambiguity, therefore multiple subband differential interferometry figure can be formed to a storehouse, then by least square method, estimate the elevation of each pixel:

$\underset{h}{\mathrm{Min}}\mathrm{\&Sigma;}{(h-\frac{{\mathrm{\&Delta;\&phi;}}_i}{2\mathrm{\&pi;}}{H}_{2\mathrm{\&pi;},i})}^2---\left(13\right)$

In above formula, △ φ i represents pixel i phase place in a certain interferogram in storehouse, H _{2 π, i}the elevation blur level that represents this point, h represents the estimated value of this point height;

B) measurement of Surface Deformation Field

Because the deformational displacement blur level of subband differential interferometry figure is very large, therefore while being converted to deformation quantity by differential phase hardly with carrying out phase unwrapping, but in order to eliminate atmospheric effect, vertical error and other various noises, still be necessary to utilize the time series analysis of one group of subband differential interferometry figure, particularly short baseline set method, tries to achieve accurate Surface Deformation Field and dynamic change thereof; It is identical substantially that the Time Series Analysis Method of subband differential interferometry figure is interfered time Sequence Analysis Method with traditional full bandwidth data, unique different be the choosing method of high coherent point; Here utilized two criterions to determine candidate's high coherent point:

$\left\{\begin{array}{c}\frac{\mathrm{\&sigma;}}{\mathrm{\&mu;}}<c\\ \mathrm{SCR}>d\end{array}\right.---\left(14\right)$

In above formula, μ and σ are average and the variances of pixel intensity level in storehouse time range, c is the threshold value that the ratio of variance and average is set, and it can weigh this pixel whether stablize of amplitude in storehouse time range substantially, and generally value is 0.2 left and right; SCR is that the signal to noise ratio of pixel is Signal to Clutter Ratio, and in subband differential interferometry figure, the signal to noise ratio of carrying out the pixel of time series analysis must be greater than a threshold value d, general value 25dB; Chosen after candidate point, utilized the speed of general period map method or short baseline set method acquisition Surface Deformation Field, linear and non-linear deformation quantity; If by artificial reverberator, the sudden change of large scale and the measurement of non-linear deformation can not be restricted.

Images