CN103728604B - A kind of broadband synthetic aperture radar sub-band interferometric data disposal route - Google Patents

Abstract

A kind of broadband synthetic aperture radar sub-band interferometric data disposal route, the method has six large steps: step one: the self-adapting precision registration of major-minor image and the resampling of sub-picture; Step 2: the optimization of subband spectrum position and bandwidth; Step 3: the bandpass filtering of major-minor image; Step 4: subband interferogram and the map generalization of subband differential interferometry; Step 5: the phase compensation of reference surface and landform altitude (DEM); Step 6: the storehouse process of subband differential interferometry figure.The present invention, in most of the cases without the need to phase unwrapping, can obtain absolute radar interference phase place, thus can directly obtain digital topography map or Ground Deformation figure by conversion, eliminates the error diffusion because phase unwrapping mistake causes; It is specially adapted to the digital terrain mapping of alpine region and the landslide of large deformation speed and cave-in areas and monitors.

Description

A kind of broadband synthetic aperture radar sub-band interferometric data disposal route

Technical field

The present invention relates to a kind of broadband synthetic aperture radar sub-band interferometric data disposal route, belong to interfering synthetic aperture radar technical field.It in most of the cases without the need to phase unwrapping, can obtain absolute radar interference phase place, thus can directly obtain digital topography map or Ground Deformation figure by conversion.

Background technology

Interfering synthetic aperture radar is one and extracts the technology of terrain object three-dimensional spatial information according to the phase data of high-resolution radar image.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 all be wound around by the angle of 2 π, so-called phase unwrapping, is exactly by the phase gradient between accumulation radar interference image primitive, thus obtains the process of pixel absolute phase.This process is wanted to obtain correct result, and for topographical surveying, a necessary condition is, in spatial domain, the absolute value of the phase differential between adjacent two pixels must not be greater than π; For the distortion measurement of Time Series Analysis Method utilizing 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 alpine region, or to have carry out interferometry compared with the landslide of large deformation velocity gradient and the Deformation Field in mine time, the condition that the absolute value of phase differential above-mentioned must not be greater than π usually can not be met, thus causing measurement result to occur very big error, sometimes even there will be the situation without separating.Theoretical according to radar interference, the problem of this phase unwrapping be solved, will try every possible means to improve interferometric phase exactly 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 π}i.e. elevation blur level, it represents the ground elevation that interference phase difference is corresponding when 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 π}i.e. deformation blur level, it represents the Ground Deformation amount that differential interferometry phase differential is corresponding when being 2 π.

Obviously, as can be seen from (1) and (2) formula, increase elevation and deformation blur level that radar carrier wavelength can increase 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 utilizes the characteristic of wideband-radar signal, in frequency domain, first radar image is split into several sub-band images, then the radar sub-band images of correspondence is carried out interference computing, height two subband interferograms 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, thus can the wavelength of original centimetre-sized be expanded 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 broadband synthetic aperture radar sub-band interferometric data disposal route, which overcome the deficiencies in the prior art, in most of the cases without the need to phase unwrapping, absolute radar interference phase place can be obtained, thus can directly obtain digital topography map or Ground Deformation figure by conversion, eliminate the error diffusion because phase unwrapping mistake causes.

2. technical scheme: the present invention is a kind of broadband synthetic aperture radar sub-band interferometric data disposal route, and the method concrete steps are as follows:

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

The measurement of higher degree or have mainly being applied in alpine region due to subband interference technique is measured compared with the landslide of large deformation speed and cave-in areas, when major-minor image carries out registration, not only to consider radar interference system (imaging geometry, interfere geometry, radar orbit etc.) skew of the position of same target on major-minor image, ground that causes, also to consider the position skew on the 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 carry out resampling, 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) illustrate in sub-picture be positioned at the i-th row jth row pixel relative to the side-play amount of pixel position same in master image, it be the side-play amount that causes because of elevation or deformation of system deviation and target itself with.Displacement when offset (i, j) gets maximal value by the covariance function of the neighborhood of pixel and the neighborhood of corresponding master image p pixel that are positioned at (i, j) in sub-picture s 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 represents major-minor image respectively, c represents in major-minor image with the covariance function of two windows of the capable M row of the N at the i-th HangjLie Wei center, n and m represents row and column side-play amount respectively, k and l represents in window the column locations participating in 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, and radar image resolution loss can be made so 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 in secondary interferogram in the form of noise.And according to signal processing theory, ideal bandpass filter is not attainable, that is, 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 away better.Obviously, two perfect conditions can not be met simultaneously.The selection of a compromise, be that the bandwidth bw of upper and lower subband equals 1/3rd of radar original signal bandwidth B, upper and lower subband is positioned at f _c± f ₀position on, namely

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

Under such optimum choice, the ratio of simulating wavelength and carrier wavelength is:

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

Generally:

$\left\{\begin{array}{c}{f}_0=\frac{B-\mathrm{bw}}{2}\\ {\mathrm{\λ}}^{\′}=\frac{f_c}{B-\mathrm{bw}}\·\mathrm{\λ}\end{array}\right.---\left(8\right)$

In three formulas above, B is radar pulse bandwidth, and fc represents carrier frequency, and λ represents carrier wavelength, and bw represents upper and lower subband bandwidth, f ₀represent the side-play amount of upper and lower subband center frequency for 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 respective complex image of subband up and down.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, namely-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 selects 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 the constant that control window leads to bandwidth sum stopband attenuation tradeoff, usually gets between 3 ~ 6.

Fig. 2, Fig. 3 and Fig. 4 respectively illustrate original radar SLC image and by the frequency spectrum on oblique distance (SlantRange) 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, generates subband differential interferometry figure.Note P and S is respectively major-minor image, and PS is their subband differential interferometry figure, then

$\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 the figure of subband up and down split off by major-minor image respectively, PS represents the subband differential interferometry figure be made up 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 SAR Interferometry, the WGS84 spheroid surface of the earth also need to superpose survey the DEM in region, 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 the distance R calculating 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 only compensate for reference face phase place, get h=0, a and b is 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, by the poor Δ R of the distance of pixel in major-minor, just can calculate the analogue phase of reference surface in subband differential interferometry figure and 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 interference.Fig. 5 is the schematic diagram to the main process of the formation of subband differential interferometry figure from the radar image after two registrations.Utilize (1) formula and (2) subband differential interferometry figure just can be converted into earth's surface elevation or Surface Deformation Field by the elevation expressed by formula and deformation blur level.

Step 6: the storehouse process of subband differential interferometry figure

A) generation of Digital height model (DEM)

Although the elevation blur level of subband differential interferometry figure has been enhanced a lot of times, also the variance of the error of interferometric phase has also been expanded simultaneously and be close to identical multiple, the phase place of individual subband differential interferometry figure can containing larger noise.Because the phase place of interferogram is proportional to elevation by the factor of phase ambiguity, therefore multiple subband differential interferometries figure can be formed a storehouse, then use the elevation of each pixel of Least Square Method:

$\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 the phase place of pixel i in storehouse in a certain interferogram, H _{2 π, i}represent the elevation mould of 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 hardly with carrying out phase unwrapping when being converted to deformation quantity by differential phase, but in order to eliminate atmospheric effect, vertical error and other various noises, still the time series analysis utilizing one group of subband differential interferometry figure (storehouse) is necessary, particularly Short baseline collection method, tries to achieve accurate Surface Deformation Field and dynamic change thereof.The Time Series Analysis Method of subband differential interferometry figure is identical with traditional full bandwidth data interferes time Sequence Analysis Method substantially, uniquely unlike the choosing method of high coherent point.In the present invention, make use of two criterions to determine the high coherent point of candidate:

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

In above formula, μ and σ is average and the variance of pixel intensity level within the scope of stack time, c is a threshold value of the ratio set to variance and average, and it can weigh this pixel whether the stablizing of amplitude within the scope of stack time substantially, and general value is about 0.2.SCR is the signal to noise ratio (SignaltoClutterRatio) of pixel, and in subband differential interferometry figure, the signal to noise ratio can carrying out the pixel of time series analysis must be greater than a threshold value d, general value 25dB.After have chosen candidate point, utilize general period map method or Short baseline collection method can obtain the speed of Surface Deformation Field, linearity and non-linearity deformation quantity.If by artificial reverberator, the sudden change of large scale (several tens cm is meter level even) and the measurement of non-linear deformation all can not be restricted in the present invention.

3. advantage and effect: the invention provides a kind of broadband synthetic aperture radar sub-band interferometric data disposal route, its advantage is:

(1) this method utilizes the subband signal of wideband radar to interfere, hundred times of ground are become to increase analog carrier wavelength and interferometry elevation blur level, in alpine region with a varied topography, the solution process of twining of mistake very easily need not be caused can to obtain radar interference absolute phase, thus obtain the digital elevation on earth's surface.

(2) this method utilizes the subband signal differential interferometry of wideband radar, directly can obtain earth's surface the several tens cm even deformation of meter level, greatly expand the measurement range of conventional differential interferometer radar.

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

(4) because this method constructs the subband signal of two different center frequency, thus as double-frequency GPS receiver, the change of ionospheric electron density can be measured, eliminates the phase error that this change causes, improve 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 subband center frequency is relative to the side-play amount of 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 interferes process flow diagram

Fig. 6. the TerraSAR-X radar image of Shu Ping slide area

Fig. 7. No. 8 and No. 9 reverberator slide displacement amounts three kinds of measuring method results contrast figure on the sliding mass of tree level ground

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

Embodiment

For landslide, tree level ground, Three Gorge Reservoir area deformation monitoring, the concrete operation step of the present invention in practical engineering application is described.See Fig. 8, the present invention is a kind of broadband synthetic aperture radar sub-band interferometric data disposal route, and the method concrete steps are as follows:

Step one: 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 to generate for the purpose of DEM, be only applicable to alpine region, for deformation monitoring, be only applicable to the area with large deformation speed.No matter which kind of object, all requires that ground monitored (or mapping) target can have higher radar coherence.For the deformation monitoring of the high coherent point by time Sequence Analysis Method, then require that ground monitored (or mapping) target can have higher radar signal to noise ratio (SCR>25dB), select if desired and artificial reverberator is installed.Greatly, long-term crops and fruit tree cover the landslide morphology change of tree level ground, Three Gorges, and landslide unstable region slippage is very large, and conventional interference radar cannot monitor out the deformation data on this landslide, therefore should adopt the artificial reverberator of installation and subband differential interference method.Tree level ground sliding mass installs 14 artificial reverberators, has separately had 4 to be arranged on outside sliding mass, as the benchmark measured and verify.

2. data decimation

Choose the spaceborne wideband synthetic aperture radar SLC image that several cover workspace, more than bandwidth requirement 100MHz.For this purpose, have chosen the 29 width SLC data on the April 9,22 days to 2010 April in 2009 of the spaceborne X-band of German TerraSAR-X 3 meters of resolution.Data fabric width 30 kilometers, 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 date are in shown in table 1 and table 2:

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: the radar data date is shown

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 reference 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, again by means of artificial reverberator in the landslide monitoring of tree level ground, therefore choosing of master image is not subject to chatted restriction above.In order to time series analysis is convenient, choose first image, namely the data on January 7th, 2012 are master image.Based on (1) formula, it is vectorial relative to the relativity shift of master image that (2) formula of utilization and (3) formula calculate all pixels in all the other each sub-pictures, with the Sinc function of displacement, self-adaption two-dimensional resampling is carried out to each sub-picture, 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)$

Section 1 in formula (1) represents in master image the departure of the position being arranged in the pixel of the i-th row jth row and the same place of sub-picture, and this departure is a bivector, i.e. Section 2 in formula (1).Bias vector is that the component that the deformation of terrain object that in the component and master image produced by system, this pixel is corresponding own produces is formed, and to be respectively in formula (1) shown in third and fourth.

When the parameters such as satellite orbit number and radar imagery time are known, the offset component that system produces first can be calculated.Based on this systematic component, roughly can 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, N are the ranks number of window, p and s represents major-minor image respectively.M and n when making covariance function obtain maximal value is the side-play amount of pixel (i, j) in sub-picture in master image, as the formula (3).

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

At radar carrier frequency f _cwhen known with 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 is at most significant end, can increase bandpasstilter stopband decay so to greatest extent, thus reduce the noise in subband differential interferometry figure.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 reduces three times, obvious resolution can increase three times.After bw determines, f ₀=(B-bw)/2, the simulating wavelength of subband differential interferometry figure is the k of radar carrier wavelength doubly, k=f _c/ (B-bw).According to landform or the deformation data maximal value of workspace, check that whether elevation blur level or the deformation blur level of subband interference be suitable, otherwise, the size of k can be adjusted by the size changing bandpass filter centre frequency and bandwidth, thus 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 the optimal design value of bandpass filter below like this:

$\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, bring the parameter in formula (4) into formula (5), get β=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 window to lead to the compromise of bandwidth sum stopband attenuation

A constant of relation, gets between 3 ~ 6 usually

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, the s after all 29 registrations ₁..., s ₂₈) convert frequency domain figure picture (P, S to by the quick Fourier transforms (2D-FFT) of two dimension ₁..., S ₂₈).Two band-pass filter up and down that these 29 frequency domain figure pictures design respectively by step 4, generate the upper sub-band images in respective frequency domain and lower sub-band images (P _up, P _low, S _{1_up}, S _{1_low}..., S _{28_up}, S _{28_low}).Before filtering, sub-band images will pass through demodulation, is also spectrum shift, and by logical for band spectrum shift in base band, the analog carrier frequency making 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 ensure the correct of the same and result of the frequency spectrum of follow-up signal process.And then utilize two dimension to convert these 58 frequency domain figure pictures to 58 spatial domain image (p against Fourier transforms (2D-IFFT) fast _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, and Fig. 2, Fig. 3 and Fig. 4 respectively illustrate original radar SLC image and by the frequency spectrum on the later oblique distance of sub-band images up and down (SlantRange) direction of band-pass filter.

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

With on April 22nd, 2009 image for master image, all the other are all the mode of sub-picture, generate 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 generation subband differential interferometry figure, and the interference image that it is produced by the lower sub-band images of right two of this interference and the interference images that two upper sub-band images produce again are interfered and produced.Because the bandwidth of now subband differential interferometry figure only has 1/3rd original (when optimal design), pixel resolution has reduced three times, thus do again interfere time, the multiple look processing of at least 3x3 should be done, to reduce speckle noise.

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

The phase compensation procedure of reference surface and landform altitude and conventional differential interfere in phase compensation 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 in face, if do Ground Deformation monitoring, except making reference except the phase compensation in face, also needs the phase compensation doing landform altitude.In order to the orbit data match with the satellite as radar carrier, choose the surface of earth WGS84 ellipsoid model as reference 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 calculate with (10) formula, interfere different from conventional differential, at this moment guinea pig carrier wavelength is k times of primary carrier wavelength, and the value of k is shown in formula (11).The simulation interferometric phase that this is produced by reference surface and landform altitude is deducted from subband differential interferometry figure.If only estimate landform altitude, the h in (9) formula is made to be 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 only compensate for reference face phase place, get h=0, a and b is 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, by the poor Δ R of the distance of pixel in major-minor, just can calculate the analogue phase of reference surface in subband differential interferometry figure and 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 interference:

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

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

Step 6: the storehouse process of subband differential interferometry figure

1. the filtering of subband differential interferometry figure

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

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

Because deformation (displacement) blur level of subband differential interferometry figure is very large, therefore hardly with carrying out phase unwrapping when being converted to deformation quantity by differential phase, but in order to eliminate atmospheric effect, vertical error and other various noises, still the time series analysis utilizing one group of subband differential interferometry figure (storehouse) is necessary, particularly Short baseline collection method, tries to achieve accurate Surface Deformation Field and dynamic change thereof.The Time Series Analysis Method of subband differential interferometry figure is identical with traditional full bandwidth data interferes time Sequence Analysis Method substantially, uniquely unlike the choosing method of high coherent point.In the present invention, make use of two criterions to determine the high coherent point of candidate:

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

In above formula, μ and σ is average and the variance of pixel intensity level within the scope of stack time, c is a threshold value of the ratio set to variance and average, and it can weigh this pixel whether the stablizing of amplitude within the scope of stack time substantially, and general value is about 0.2.SCR is the signal to noise ratio (SignaltoClutterRatio) of pixel, and in subband differential interferometry figure, the signal to noise ratio can carrying out the pixel of time series analysis must be greater than a threshold value d, general value 25dB.After have chosen candidate point, utilize general period map method or Short baseline collection method can obtain the speed of Surface Deformation Field, linearity and non-linearity deformation quantity.If by artificial reverberator, the sudden change of large scale (several tens cm is meter level even) and the measurement of non-linear deformation all 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, read these 18 artificial reverberators and often opening the phase place in subband differential interferometry figure, and deduct the phase place of a reference point.In landslide, tree level ground, No. 6 artificial reverberators are chosen to be reference point.After obtaining the phase place of each artificial reverberator in each subband differential interferometry figure, by formula (13), the displacement on the radar ray direction of each artificial reverberator within the corresponding interferogram time period can be calculated:

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

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

3. the subsequent treatment of subband differential interferometry figure

The solution of subband differential interferometry figure twines, phase place-elevation conversion, 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 this few part is not within the scope of the present invention.What make an exception is only the carrier wavelength in (14) and (15) formula blur level expression formula is made 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 π}i.e. elevation blur level, it represents the ground elevation that interference phase difference is corresponding when 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 π}i.e. deformation blur level, it represents the Ground Deformation amount that differential interferometry phase differential is corresponding when being 2 π.

Claims (1)

1. a broadband synthetic aperture radar sub-band interferometric data disposal route, is characterized in that: the method concrete steps are as follows:

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

The measurement of higher degree or have mainly being applied in alpine region due to subband interference technique is measured compared with the landslide of large deformation speed and cave-in areas, when major-minor image carries out registration, not only to consider the skew of the position of same target on major-minor image, ground that radar interference system causes, also will consider the position skew on the major-minor image that terrain object causes because of larger elevation or larger deformational displacement; Here, sub-picture utilizes Sinc function to carry out resampling, and the displacement of Sinc function is determined by following formula:

$Offset(i,j)=(\mathrm{\&Delta;}i,\mathrm{\&Delta;}j)={\left({\left(\mathrm{\&Delta;}i,\mathrm{\&Delta;}j\right)}_{i,j}\right)}_{System}+{\left({\left(\mathrm{\&Delta;}i,\mathrm{\&Delta;}j\right)}_{i,j}\right)}_{Object}---\left(3\right)$

In above formula, offset (i, j) illustrate in sub-picture be positioned at the i-th row jth row pixel relative to the side-play amount of pixel position same in master image, it be the side-play amount that causes because of elevation or deformation of system deviation and target itself with; Displacement when offset (i, j) gets maximal value by the covariance function of the neighborhood of pixel and the neighborhood of corresponding master image p pixel that are positioned at (i, j) in sub-picture s determines, sees (4) and (5) formula below

$c_{ps}{(n,m)}_{i,j}=\underset{k=-N/2}{\overset{k=N/2}{\&Sigma;}}\underset{l=-M/2}{\overset{l=M/2}{\&Sigma;}}\left(p\right(k,l)-\tilde{p})\&CenterDot;\left(s\right(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 represents major-minor image respectively, c represents in major-minor image with the covariance function of two windows of the capable M row of the N at the i-th HangjLie Wei center, n and m represents row and column side-play amount respectively, k and l represents in window the column locations participating in correlation computations pixel, (Δ i, Δ j) represents the ranks number when c obtains maximal value, i.e. (n _max, m _max);

Step 2: the optimization of subband spectrum position and bandwidth

In ideal conditions, wish that subband bandwidth is wide as much as possible, radar image resolution loss can be made so 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 in secondary interferogram in the form of noise; And according to signal processing theory, ideal bandpass filter is not attainable, that is, 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 away better; Obviously, two perfect conditions can not be met simultaneously, and the selection of a compromise, be that the bandwidth bw of upper and lower subband equals 1/3rd of radar original signal bandwidth B, upper and lower subband is positioned at f _c± f ₀position on, namely

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

Under such optimum choice, the ratio of simulating wavelength and carrier wavelength is:

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

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

In three formulas above, B is radar pulse bandwidth, f _crepresent carrier frequency, λ represents carrier wavelength, and bw represents upper and lower subband bandwidth, f ₀represent the side-play amount of upper and lower subband center frequency for 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 respective complex image of subband up and down; 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, namely-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 selects Caesar's window and 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 the constant that control window leads to bandwidth sum stopband attenuation tradeoff, 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, generates subband differential interferometry figure; Note P and S is respectively major-minor image, and PS is their subband differential interferometry figure, then

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

In above formula, P _up, P _low, S _up, and S _lowrepresent the figure of subband up and down split off by major-minor image respectively, PS represents the subband differential interferometry figure be made up 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 SAR Interferometry, the WGS84 spheroid surface of the earth also need to superpose survey the DEM in region, 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 the distance R calculating 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 only compensate for reference face phase place, get h=0, a and b is 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, by the poor Δ R of the distance of pixel in major-minor, just calculate the analogue phase of reference surface in subband differential interferometry figure and 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 interference, utilizes formula and subband differential interferometry figure is just converted into earth's surface elevation or Surface Deformation Field by the elevation expressed by formula and deformation blur level;

Step 6: the storehouse process of subband differential interferometry figure

A) generation of Digital height model and DEM

Although the elevation blur level of subband differential interferometry figure has been enhanced a lot of times, also the variance of the error of interferometric phase has also been expanded simultaneously and be close to identical multiple, the phase place of individual subband differential interferometry figure can containing larger noise; Because the phase place of interferogram is proportional to elevation by the factor of phase ambiguity, therefore multiple subband differential interferometries figure is formed a storehouse, then uses the elevation of each pixel of Least Square Method:

$\underset{h}{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 the phase place of pixel i in storehouse in a certain interferogram, H _{2 π, i}represent the elevation blur level of 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 need not phase unwrapping be carried out when being converted to deformation quantity by differential phase, but in order to eliminate atmospheric effect, vertical error and other various noises, also to utilize the time series analysis of one group of subband differential interferometry figure, try to achieve accurate Surface Deformation Field and dynamic change thereof; The Time Series Analysis Method of subband differential interferometry figure is identical with traditional full bandwidth data interferes time Sequence Analysis Method, uniquely unlike the choosing method of high coherent point; Here make use of two criterions to determine the high coherent point of candidate:

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

In above formula, μ and σ is average and the variance of pixel intensity level within the scope of stack time, and c is a threshold value of the ratio set to variance and average, and it can weigh this pixel whether the stablizing of amplitude within the scope of stack time, and value is 0.2; SCR is signal to noise ratio and the SignaltoClutterRatio of pixel, 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, value 25dB; After have chosen candidate point, period map method or Short baseline collection method is utilized to obtain the speed of Surface Deformation Field, linearity and non-linearity deformation quantity; If by artificial reverberator, the sudden change of large scale and the measurement of non-linear deformation all can not be restricted.