CN104123464A - Method for inversion of ground feature high elevation and number of land subsidence through high resolution InSAR timing sequence analysis - Google Patents

Abstract

A method for inversion of the ground feature high elevation and the number of land subsidence through high resolution InSAR timing sequence analysis includes the first step of high resolution SAR data selection, the second step of initial elevation phase estimation, the third step of selection of InSAR timing sequences, analysis of an interference image pair data set and extraction of coherent target candidate points, the fourth step of elevation phase correction and linear deformation component estimation, and the fifth step of deformation sequence estimation of coherent targets. According to the method, the defect that no high resolution DEM is available for removal of high-rise building additional phases in urban area high resolution InSAR deformation monitoring data processing can be overcome, the additional phases and deformation phases of urban high-rise buildings in the high resolution InSAR are effectively separated, and therefore monitoring precision of the high resolution InSAR in urban area land subsidence is greatly improved.

Description

A kind of method of high resolving power InSAR time series analysis inverting atural object elevation and ground settlement

Technical field

The present invention relates to a kind of method of high resolving power InSAR time series analysis inverting atural object elevation and ground settlement, belong to interfering synthetic aperture radar measuring technique (InSAR) field.It can by urban skyscraper thing, the additive phase in high resolving power InSAR interferogram separates with deformation phase place effectively, can greatly improve the monitoring accuracy of high resolving power InSAR in city ground sedimentation.

Background technology

Utilize InSAR to solve elevation (DEM) or deformation quantity, phase unwrapping is one of necessary step, this step is called solution and twines in differential interferometry is processed, in Time Series Analysis Method, be called parameter estimation, its essence is and solve the complete cycle number that is wound around phase place, basic process is the phase differential solving between consecutive point, then carries out integration according to specific path or network with certain constraint condition, the solution that solves observation scope entirety twines phase place, and then inverting Deformation Field or elevation.The prerequisite of phase unwrapping is that interferogram continuous distribution and variation are mild, meets the constraint condition that phase differential is less than π, and the whole field of behaviour is irrotational field, and fruit does not become with path disentanglement.And in fact interferogram is subject to the impact of noise or discontinuous phase place (as discontinuous abrupt slope), is difficult to meet solution to twine condition, thereby connects after need to solving according to given path again.Under high resolving power condition, thing caused phase place variation in urban skyscraper is similar to discrete abrupt slope phase place.For solving elevation, fringe density is intensive with the increase of baseline, and existence mix with deformation striped may.Thereby, how to resolve elevation and deformation phase place simultaneously, the precision that improves deformation quantity estimation is that high resolving power InSAR applies the main bugbear facing.

Landform phase compensation is the basic step that in InSAR time series analysis processing procedure, differential phase is asked for, and side by side object height journey is also the main calculation result of InSAR time series analysis.The pixel that shows as one or several point targets in intermediate resolution SAR data shows with several point targets in high-resolution radar data, the scattering center difference of each scatterer, as a solitary building, in show as a point target in point SAR image, and show with multiple different scatterers in high score SAR image.Direct method to the phase compensation of high resolving power InSAR landform is to obtain the high-precision dem data of each scattering center, utilizes dem data that Lidar, TanDEM-X etc. obtain to realize the elevation phase compensation of approximate equal resolution differential interferometry processing.But in a lot of situations of this high-resolution dem data, be difficult to obtain.Therefore, this method that compensates whole pixel elevation by outside DEM (as SRTM) will be difficult to be suitable in high-resolution data is processed, and need to study the method for single target elevation accurate Calculation.

In InSAR interferogram, as shown in Figure 1, under equal elevation, baseline is shorter for the frequency of interference fringe and the relation of baseline, and striped is more sparse, is beneficial to phase unwrapping or parameter estimation.When local deformation turns to 2 π, corresponding difference of elevation is:

${\mathrm{\Δh}}_{2\mathrm{\π}}=\frac{\mathrm{\λ}}{2}\frac{R\sin \mathrm{\θ}}{B_{\⊥}}---\left(1\right)$

Here Δ h, _{2 π}for elevation blur level, characterize the sensitivity that interferometric phase changes landform altitude.The impact changing with respect to elevation, deformation phase place and elevation are irrelevant, but are subject to the impact that elevation changes.

According to the relation of deformation phase place and elevation phase place, known elevation changes relevant to baseline on the impact of deformation phase place:

$\mathrm{\δd}=\frac{B_{\⊥}\·\mathrm{\δh}}{R\sin \mathrm{\θ}}---\left(2\right)$

With Cosmo-skymed data instance, Fig. 2 has characterized the impact of vertical error on deformation quantity.In the time that vertical parallax is 50m, the phase place variation that the difference of elevation of 50m causes is less than 0.41 complete cycle, is 0.82 complete cycle when 100m.Under same condition, when difference of elevation is 100m, the caused striped of 100m baseline is changed to 1.64 complete cycles, is while being 50m 2 times of baseline.Obviously, under short base line condition, fringe density reduces, and this is conducive to the accurate extraction to deformation data.In addition, be conducive to elevation inversion accuracy with respect to long baseline, the present invention considers to fully utilize length base line interference figure iteration correction elevation estimated value, to reach the object that improves deformation estimated accuracy.

Summary of the invention

1. object: a kind of method that the object of this invention is to provide high resolving power InSAR time series analysis inverting atural object elevation and ground settlement.It can overcome in the high resolving power InSAR Deformation Monitoring Data processing of city removes high-rise additive phase without high resolution DEM, by urban skyscraper thing, the additive phase in high resolving power InSAR separates with deformation phase place effectively, thereby can greatly improve the monitoring accuracy of high resolving power InSAR in city ground sedimentation.

2. technical scheme: the present invention is a kind of method of high resolving power InSAR time series analysis inverting atural object elevation and ground settlement, and the method concrete steps are as follows:

Step 1: high resolution SAR data decimation

High-resolution radar satellite system taking TerraSAR-X and COSMO-Skymed as representative provides data source for carrying out InSAR inverting city atural object elevation and the deformation monitoring that becomes more meticulous.High resolving power InSAR is with respect to intermediate resolution InSAR technology, and its overall advantage is embodied in two aspects, i.e. (I) high density coherent point target and short period (4-16 days); (II) accurate location to ground point target.The data decimation of high resolving power InSAR data processing wants to cover as far as possible on meeting spatial the orbital data in complete city, and data can receive continuously in time.

Step 2: initial elevation phase estimation

In the constructed two-dimensional parameter estimation model of InSAR deformation time series analysis, consider the spatial coherence of atmosphere, adjacent 2 points (being less than atmosphere correlation distance) are asked to poor to weaken the impact of atmosphere.The mutual deviation of Coherent Targets i and j differential interferometry phase place is:

$\mathrm{\Δ}{\mathrm{\φ}}_{i,j}^k=[C_B\·B^{\left(k\right)}\·{\mathrm{\Δ\ϵ}}_{i,j}+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·T^{\left(k\right)}\·{\mathrm{\Δv}}_{i,j}]+{\mathrm{\μ}}_{\mathrm{NL}}^{\left(k\right)}+{\mathrm{\α}}^{\left(k\right)}+n^{\left(k\right)}---\left(3\right)$

In above formula, C _bfor the coefficient relevant to vertical parallax, T is time basis, and Δ ε is relative altitude error, and Δ v is relative deformation speed, μ _nLfor non-linear deformation quantity, α is atmospheric phase, and n is noise, and k represents interferogram number (relevant with the combination of interferogram sequence).Establishing target function is as follows:

${\mathrm{\φ}}_{\mathrm{mod}\mathrm{el}}(i,j{,T}^{\left(k\right)})=C_B{\·B}^{\left(k\right)}\·{\mathrm{\Δ\ϵ}}_{i,j}+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·T^{\left(k\right)}\·{\mathrm{\Δv}}_{i,j}---\left(4\right)$

Above formula is deducted from phase place mutual deviation formula (3), obtains residual phase and be:

${\mathrm{\Δw}}_{i,j}^k={\mathrm{\Δ\φ}}_{i,j}^{\left(k\right)}-[C_B\·B^{\left(k\right)}\·{\mathrm{\Δ\ϵ}}_{i,j}+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·T^{\left(k\right)}\·{\mathrm{\Δv}}_{i,j}]={\mathrm{\μ}}_{\mathrm{NL}}^{\left(k\right)}+{\mathrm{\α}}^{\left(k\right)}+n^{\left(k\right)---\left(5\right)}$

Obviously,, when the parameter Δ ε of objective function and Δ v are in the time accurately estimating, residual phase will minimize.

The present invention takes the estimation of being undertaken on the basis of differential interferometry figure phase unwrapping, this up-to-date style (3) is converted to two-dimensional linear function, can be by setting up the Delanay triangulation network or utilizing redundant network to build more complicated annexation strengthening and treat the constraint of resolving system of equations, the Coherent Targets that utilizes contiguous rule that all distances are met to atmosphere correlation distance couples together, after the mutual deviation having solved between consecutive point, solve height value and the rate of deformation field of each target with respect to reference point by least square or average weighted method.

When enough hour of the interferogram time interval and Space Baseline, distortion phase place can be ignored, and two dimensional model (formula (3)) can be transferred to one-dimensional model (formula (6)), and then resolved and obtain DEM.

$\mathrm{\Δ}{\mathrm{\φ}}_{i,j}^k=C_B\·B^{\left(k\right)}\·{\mathrm{\Δ\ϵ}}_{i,j}+{\mathrm{\α}}^{\left(k\right)}+n^{\left(k\right)}---\left(6\right)$

(1) choose ultra-short Time and Space Baseline interference image to data set

According to above-mentioned theory, according to the relation of the height of buildings, interferogram baseline and landform phase place fringe density, choose the average depth of building in study area and retrain as elevation, determine that the interference image needing while calculating object height journey is primitively to sequence.Access time baseline and Space Baseline all shorter interference image to sequence as the data set initially solving, set initial time and Space Baseline threshold value and be respectively T ₁and B ₁to the interference image of data centralization to interfering processing, carry out landform phase correction without outside DEM, in the present invention, fully utilize the method that Fast Fourier Transform (FFT) (FFT) is estimated and fitting of a polynomial is estimated and remove orbit error and the tendency interference fringe in interferogram.Generate ultrashort space-time base line interference diagram data collection, coherence map data set and intensity map data set.Each interferogram is carried out to phase unwrapping, according to one-dimensional model repeatedly iteration correction vertical error phase estimation obtain original DEM.It should be noted that, operations all before geocoding are all carried out under radar fix system.

(2) utilize point target recognition methods to extract Coherent Targets candidate point

For the interferogram sequence of above-mentioned all generations, in the present invention, comprehensive employing amplitude dispersion index (Amplitude Dispersion Index) and coefficient of coherence (coherence) screen Coherent Targets candidate point.

The computing formula of amplitude dispersion index is:

$D_A=\frac{{\mathrm{\σ}}_A}{m_A}---\left(7\right)$

Wherein, σ _aand m _abe respectively standard deviation and the average of pixel amplitudes value.A given appropriate threshold value d _abe defined as Coherent Targets candidate point lower than the pixel of threshold value.

The Coherence Estimation formula of radar interference phase diagram is:

$\widetilde{\mathrm{\γ}}=\left|\frac{\frac{1}{N}{\mathrm{\Σ}}_{i=0}^N{M}_i{S}_i^{*}}{\frac{1}{N}{\mathrm{\Σ}}_{i=0}^N{M}_i{M}_i^{*}\frac{1}{N}{\mathrm{\Σ}}_{i=0}^N{S}_i{S}_i^{*}}\right|---\left(8\right)$

Coefficient of coherence sequence γ according to each pixel point in coherence map _iwith given coefficient of coherence threshold value if mean so this pixel is defined as to Coherent Targets candidate point.

(3) iteration correction vertical error phase estimation atural object elevation repeatedly

For each Coherent Targets point, utilize ultrashort space-time base line interference diagram data collection and according to repeatedly iteration correction vertical error phase place of one-dimensional model formula (6) Suo Shu, the vertical error correction phase place of repeatedly estimating is added, obtains the elevation at Coherent Targets point place.Then, point vector is converted to raster data, generates the original DEM of study area,

Step 3: choose InSAR time series analysis interference image to data set and extract Coherent Targets candidate point

All interference images that Space Baseline and time basis are no more than to a certain setting threshold are to interfering processing, and the original elevation artificially generated terrain phase place that step 2 is generated, for the elevation phase compensation of all interferograms, generates differential interferometry graphic sequence.For the orbit error occurring in single differential interferometry figure and tendency interference fringe, the method that comprehensive utilization Fast Fourier Transform (FFT) (FFT) is estimated and fitting of a polynomial is estimated is removed, final differential interferometry diagram data collection and corresponding coherence map data set and the intensity map data set generating for InSAR time series analysis.To the differential interferometry graphic sequence of above-mentioned all generations, comprehensively adopt amplitude dispersion index and coefficient of coherence to screen Coherent Targets candidate point, reduce Coherent Targets quantity redundancy under high resolution SAR condition.

Step 4: elevation phase correction and linear deformation component are estimated

Solution twines each the width differential interferometry figure in InSAR time series analysis differential interferometry graphic sequence.Access time and Space Baseline threshold value are respectively T ₂and B ₂interference image to data set, re-use one-dimensional model secondary and estimate elevation residual value dh _ori.By above-mentioned one-dimensional model gained elevation residual value dh _orias the initial elevation of deformation time series analysis two-dimensional parameter estimation model, the step-length using Δ T, Δ B and Δ dh as time, Space Baseline and elevation correction, by T ₂+ Δ T and B ₂interferogram in+Δ B participates in the sequence of two-dimensional parameter estimation, sets initial maximum elevation correction threshold DH _{max_ori}, and estimate for the first time elevation modified value and rate of deformation.Progressively increase the interferogram quantity for iterative estimation with step delta T and Δ B respectively, and progressively reduce maximum elevation correction threshold with step delta dh, progressively revise elevation estimated value and rate of deformation by iterative processing repeatedly.All participate in calculating until meet all differential interferometry figure of time and Space Baseline and maximum elevation correction threshold, estimate final elevation correction and rate of deformation.Be the estimation to removing height value because successive iteration solves the elevation correction obtaining, thereby final atural object elevation estimated value is successive iteration modified value sum, that is:

$H=\underset{i=1}{\overset{u}{\mathrm{\Σ}}}{\mathrm{\ϵ}}^i---\left(9\right)$

Atural object elevation estimated value when H is iteration termination in formula, u is total iterations.

Step 5: the deformation sequence estimation of Coherent Targets

For thering is remarkable non-linear deformation process, still need residual phase to carry out more complicated processing, to extract non-linear deformation quantity.The prerequisite of processing is still the spatio-temporal frequency characteristic based on out of phase.The residual phase of removing from differential phase after elevation and linear velocity is:

${\mathrm{\Δw}}_{i,j}^k={\mathrm{\μ}}_{i,j}^k+a_{i,j}^k+n_{i,j}^k---\left(10\right)$

Owing to calculating complete cycle phase place, thereby residual phase is phase place main value, and its size meets:

$\left|{\mathrm{\&Delta;w}}_{i,j}^k\right|<\mathrm{\&pi;}---\left(11\right)$

Here first to realize the filtering processing to residual phase in single differential interferometry figure.Because atmosphere shows space low frequency characteristic, and the spatial variations scope of nonlinear deformation is less, and atmospheric phase shows as high pass characteristic relatively.Thereby, residual phase bitmap is carried out low pass spatial filtering and can further be weakened the impact of atmosphere.The non-linear deformation quantity solving under short baseline set condition, is different from the processing policy of PSInSAR, directly atmosphere is not estimated, but is weakened the impact of atmosphere by filtering.Residual phase is after the high-pass filtering of space, and atmosphere component weakens.Now, further processing is to solve corresponding with radar data not nonlinear deformation amount in the same time, comprises the impacts such as noise.According to the relation of major-minor image, the residual phase after phase unwrapping can be decomposed into:

${\mathrm{\&Delta;w}}_i^k={\mathrm{\&Delta;w}}_i^p-{\mathrm{\&Delta;w}}_i^q,\&ForAll;k=1,\&CenterDot;\&CenterDot;\&CenterDot;,N---\left(12\right)$

In formula, p and q represent to generate k and open the acquisition time of the major-minor image of differential interferometry figure.Owing to adopting short baseline principle, solving each moment when corresponding residual phase, there will be rank defect system of equations, the way addressing this problem is exactly singular value decomposition method (SVD).Thus, on the basis of svd, M is opened to residual phase and carries out time domain low-pass filtering treatment, to extract final non-linear deformation quantity, can be expressed as:

${\left(d_i^p\right)}_{\mathrm{NL}\_\mathrm{est}}=\frac{\mathrm{\&lambda;}}{4\mathrm{\&pi;}}\&CenterDot;{\left\{\right[{\mathrm{\&Delta;w}}_i^p{]}_{\mathrm{HP}\_\mathrm{space}}\}}_{\mathrm{TP}\_\mathrm{time}}---\left(13\right)$

Solving after non-linear deformation quantity, the deformation sequence that each Coherent Targets is corresponding is:

$d^p=v_{\mathrm{est}}\&CenterDot;(t^p-t^1)+{\left(d_i^p\right)}_{\mathrm{NL}\_\mathrm{est}},k=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,M---\left(14\right)$

Wherein, v _estobtain Linear deformation rate for resolving.

3. advantage and effect:

A kind of high resolving power InSAR time series analysis inverting atural object elevation that the present invention proposes is the same with the principle of the method for ground settlement and the InSAR deformation monitoring in other field.But feature of the present invention is, it participates in without outside DEM, directly utilizes two dimensional model Simultaneous Inversion to obtain atural object elevation and rate of deformation, thereby can greatly improve the precision of high resolving power InSAR in city ground settlement monitoring.

Brief description of the drawings

Fig. 1. the relation of fringe density and vertical parallax.For the same discrepancy in elevation, vertical parallax is longer, and striped is more intensive.

Fig. 2. with Cosmo-skymed3m data instance, the impact of vertical error on deformation quantity is discussed.Under short base line condition, fringe density reduces, and is conducive to the accurate extraction of deformation data.

Fig. 3. the process flow diagram of high resolving power InSAR time series analysis inverting atural object elevation and ground settlement.

Fig. 4. participate in the interferogram formation axis distribution plan of InSAR time series analysis.Wherein, vertical parallax is 300m, and the time interval is 1 year.

Fig. 5 (a) estimates for one-dimensional model the initial DEM obtaining, and image is the lower 90-degree rotation of radar fix system.

Fig. 5 (b) estimates the elevation residual error obtaining 2 times for one-dimensional model, and image is the lower 90-degree rotation of radar fix system.

The elevation modified value that Fig. 5 (c) obtains for the first iterative estimate of two dimensional model, image is the lower 90-degree rotation of radar fix system.

The elevation modified value that Fig. 5 (d) obtains for 2 iterative estimates of two dimensional model, image is the lower 90-degree rotation of radar fix system.

The elevation modified value that Fig. 5 (e) obtains for 5 iterative estimates of two dimensional model, image is the lower 90-degree rotation of radar fix system.

The elevation modified value that Fig. 5 (f) obtains for 9 iterative estimates of two dimensional model, image is the lower 90-degree rotation of radar fix system.

Fig. 6. the atural object elevation estimated value that iteration obtains, image is the lower 90-degree rotation of radar fix system.One-dimensional model elevation residual error and two dimensional model repeatedly the elevation modified value addition of grey iterative generation obtain final atural object elevation estimated value, and maximal value exceedes 300m.Wherein, black surround position is skyscraper position, Binjiang garden.

Fig. 7. Binjiang garden skyscraper Height Estimation value, image is the lower 90-degree rotation of radar fix system.

Fig. 8. Lujiazui, Shanghai City regional land subsidence rate diagram, the time interval: 2008/12/10-2010/06/23.

Fig. 9. near certain Coherent Targets deformation sequence Yuyuan Garden, the time interval: 2008/12/10-2010/06/23.

Embodiment

Taking Cosmo-skymed InSAR monitoring Lujiazui, Shanghai City regional land subsidence as example, the concrete operation step of the present invention in Practical Project is answered is described.As shown in Figure 3, the present invention is a kind of method of high resolving power InSAR time series analysis inverting atural object elevation and ground settling amount, and the method concrete steps are as follows:

Step 1: high resolution SAR data decimation

Choose the Cosmo-skymed3m resolution radar data that cover main city, Shanghai City, its acquisition time is year June in Dec, 2008 to 2010, and time span is a year and a half, and data set quantity is 35 scapes.Cosmo-skymed system is that Italian Ministry of National Defence and space office subsidize jointly, and Alenia space flight company is responsible for the constellation being made up of 4 radar satellites of development.Its 3 meters of resolution data fabric width 40km × 40km, orientation is to being respectively 3100MHz and 73.5MHz with distance to bandwidth, and heavily the visit cycle is 16 (4) days to satellite.Radar general data parameter is as shown in table 1, and the scanning date is referred to table 2.

Table 1: select Cosmo-skymed radar major parameter table

Carrier frequency	9.6000000e+009Hz
		Carrier wavelength	3.1cm
Pulse bandwidth	9.6000000e+007Hz
		Incident angle	40°
Distance is to Pixel size	1.249135m
		Orientation is to Pixel size	2.246403m
Distance is to sampling rate	1.2000000e+008Hz
		Radar scanning pattern	Band pattern
Data type	Haplopia plural number (SLC)

Table 2: radar data date table

20081210	20090111	20090212	20090228	20090316	20090401	20090409
							20090417	20090604	20090612	20090714	20090807	20090815	20090924
20091002	20091010	20091018	20091026	20091103	20091205	20091213
							20091221	20100122	20100207	20100223	20100311	20100319	20100327
20100404	20100412	20100428	20100506	20100514	20100615	20100623

Step 2: ultrashort space-time base line interference picture is to choosing and initial elevation phase estimation

Consider the relation of height, interferogram baseline and the landform phase place fringe density of buildings, get the average depth of building in study area and retrain as elevation, the interferogram sequence needing while determining first calculating atural object elevation.Here, according to the densely distributed feature of region, Lujiazui, Shanghai City high-rise, determine that buildings average height is 100m, when known baseline is 50m, its fringe density is 0.82 complete cycle, is 1.64 complete cycles when 100m.The accuracy of estimating to improve first elevation for reducing the error of phase unwrapping, selects the interferogram of striped in 1 complete cycle to participate in computing.Determine that time basis is less than 50 days, Space Baseline be less than in 50m totally 34 interference images to sequence as the data set initially solving, each interference image is processed and phase unwrapping interfering, and is generated corresponding coherence map and intensity map.For the orbit error occurring in each interferogram and tendency interference fringe, the method that comprehensive utilization Fast Fourier Transform (FFT) (FFT) is estimated and fitting of a polynomial is estimated is removed.Setting amplitude dispersion index and coefficient of coherence threshold value are respectively 1.5 and 0.675, and screening obtains Coherent Targets candidate point.Carry out terrain parameter estimation according to one-dimensional model, extract atural object elevation as shown in Figure 5 a.

Step 3: InSAR time series analysis interference image extracts choosing with Coherent Targets candidate point

Build interference image to sequence according to short baseline thought, by the time interval within 1 year, Space Baseline is less than the interference image of 300m to carrying out differential interferometry processing, be illustrated in figure 4 the interferogram baseline profile that participates in calculating, have 187 interference images pair, the initial DEM artificially generated terrain phase place of utilizing step 2 to obtain, removes the landform phase place in interferogram, and each differential interferometry figure is carried out to phase unwrapping.For the orbit error occurring in each interferogram and tendency interference fringe, the method that comprehensive utilization Fast Fourier Transform (FFT) (FFT) is estimated and fitting of a polynomial is estimated is removed, and finally obtains initial differential interferometry phase diagram and coherence map for time series analysis.On this basis, utilize amplitude dispersion index (threshold value 1.5) and coefficient of coherence (threshold value 0.675) to screen and obtain Coherent Targets candidate point respectively, and merging obtain final candidate point.

Step 4: elevation phase correction and linear deformation component are estimated

In access time interval 100 days, vertical parallax is 100m solution twines phase diagram, re-use one-dimensional model and solve elevation residual value.One-dimensional model is solved to the initial elevation that elevation residual value is estimated as two dimensional model, substep successive iteration completes the solution procedure of elevation and deformation data: (1) setting space baseline threshold is 100m, initial maximum elevation correction threshold is 45m, time basis is no more than to 150 days interferograms and has participated in original two-dimensional parameter estimation, its result as shown in Figure 5 c; (2) within 100 days, progressively increase interferogram quantity until time basis, as 350 days, and simultaneously gradually reduces maximum elevation correction threshold with elevation correction step-length 5m taking time basis step-length, iteration completes two-dimensional parameter and estimates; (3) on the basis completing in above-mentioned steps, progressively increase Space Baseline length until reach 300m with Space Baseline step-length 100m, repeat (2) and (3).Vertical parallax is less than 300m the most at last, and 1 year of the time interval participates in estimating final elevation correction and rate of deformation with interior interferogram (187).The elevation residual value that one-dimensional model is solved and the elevation modified value that successively solves gained are added, and try to achieve the final atural object elevation estimated value in main city, Shanghai City, Figure 6 shows that high stored building group district, the Lujiazui elevation that the present invention estimates.

In order to check the accuracy of depth of building estimated value, utilize the three-dimensional outdoor scene measurement data of Google Earth to compare building is high.Choose the high stored building group distributing in blocks shown in Fig. 6 black surround as the checked object of elevation estimated result, this building area is corresponding to the high-rise building (Fig. 7) of 7 dense distribution in the luxuriant Binjiang of generation garden, its average height reaches 150-170m, and number of floor levels is 45-50 layer, and interlamellar spacing is 3-4m.In these groups of building, many elevation buildings distribute comparatively uniformly, and be connected with other low buildings, the estimation of its height has the meaning of representative in the sequential regretional analysis of phase place, that is: be connected with other buildingss on the whole, be used for estimating elevation, simultaneously separate again between separately, representative in Height Estimation result relatively.In the three-dimensional outdoor scene of selection and Google, building clear height and height value (relative accuracy is better than 1m) compare, and extract respectively roof, the bottom of the building height value of 7 buildingss from Google earth and elevation estimation, calculate its clear height.The actual elevation of skyscraper that table 3 is 7 dense distribution in the luxuriant Binjiang of generation garden and InSAR estimate the comparison that elevation does, and both elevation mutual deviation basic controlling, in 10m, illustrate that atural object elevation estimated value of the present invention is comparatively reliable.

Step 5: land subsidence speed and deformation sequence estimation

As previously mentioned, after completing first elevation and solving, return the progressively improvement of carrying out elevation and rate of deformation by two-dimensional iteration, obtain final rate of deformation result, be converted into subsidence rate and carry out geocoding to relation according to deformation and sight line on this basis, Fig. 8 is Lujiazui, Shanghai City regional land subsidence rate diagram (2008/12/10-2010/06/23), and within the largest year, subsidence rate reaches 57.5mm.The residual phase of removing after elevation and linear velocity is carried out to time domain and spatial domain filtering processing, extract final non-linear deformation quantity.Linear and non-linear deformation quantity is added to the deformation sequence that obtains each Coherent Targets, and Fig. 9 is near the deformation sequence of certain Coherent Targets Yuyuan Garden.Fig. 1 is the schematic diagram that is related to of line density and vertical parallax; Fig. 2 is with Cosmo-skymed3m data instance, and the impact of vertical error on deformation quantity is discussed.

The end, table 3 buildings top point Differential altitude comparison

Claims (1)

1. a method for high resolving power InSAR time series analysis inverting atural object elevation and ground settlement, is characterized in that: the method concrete steps are as follows:

Step 1: high resolution SAR data decimation

High-resolution radar satellite system taking TerraSAR-X and COSMO-Skymed as representative provides data source for carrying out InSAR inverting city atural object elevation and the deformation monitoring that becomes more meticulous, high resolving power InSAR is with respect to intermediate resolution InSAR technology, its overall advantage is embodied in two aspects, i.e. (I) high density coherent point target and short period 4-16 days; (II) accurate location to ground point target; The data decimation of high resolving power InSAR data processing wants to cover as far as possible on meeting spatial the orbital data in complete city, and data can receive continuously in time;

Step 2: initial elevation phase estimation

In the constructed two-dimensional parameter estimation model of InSAR deformation time series analysis, consider the spatial coherence of atmosphere, ask poor to weaken the impact of atmosphere to adjacent 2; The mutual deviation of Coherent Targets i and j differential interferometry phase place is:

$\mathrm{\&Delta;}{\mathrm{\&phi;}}_{i,j}^k=[C_B\&CenterDot;B^{\left(k\right)}\&CenterDot;{\mathrm{\&Delta;\&epsiv;}}_{i,j}+\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\&CenterDot;T^{\left(k\right)}\&CenterDot;{\mathrm{\&Delta;v}}_{i,j}]+{\mathrm{\&mu;}}_{\mathrm{NL}}^{\left(k\right)}+{\mathrm{\&alpha;}}^{\left(k\right)}+n^{\left(k\right)}---\left(3\right)$

In above formula, C _bfor the coefficient relevant to vertical parallax, T is time basis, and Δ ε is relative altitude error, and Δ v is relative deformation speed, μ _nLfor non-linear deformation quantity, α is atmospheric phase, and n is noise, and k represents interferogram number, relevant with the combination of interferogram sequence; Establishing target function is as follows:

${\mathrm{\&phi;}}_{\mathrm{mod}\mathrm{el}}(i,j{,T}^{\left(k\right)})=C_B{\&CenterDot;B}^{\left(k\right)}\&CenterDot;{\mathrm{\&Delta;\&epsiv;}}_{i,j}+\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\&CenterDot;T^{\left(k\right)}\&CenterDot;{\mathrm{\&Delta;v}}_{i,j}---\left(4\right)$

Above formula is deducted from phase place mutual deviation formula (3), obtains residual phase and be:

${\mathrm{\&Delta;w}}_{i,j}^k={\mathrm{\&Delta;\&phi;}}_{i,j}^{\left(k\right)}-[C_B\&CenterDot;B^{\left(k\right)}\&CenterDot;{\mathrm{\&Delta;\&epsiv;}}_{i,j}+\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\&CenterDot;T^{\left(k\right)}\&CenterDot;{\mathrm{\&Delta;v}}_{i,j}]={\mathrm{\&mu;}}_{\mathrm{NL}}^{\left(k\right)}+{\mathrm{\&alpha;}}^{\left(k\right)}+n^{\left(k\right)---\left(5\right)}$

Obviously,, when the parameter Δ ε of objective function and Δ v are in the time accurately estimating, residual phase will minimize;

Take the estimation of being undertaken on the basis of differential interferometry figure phase unwrapping, this up-to-date style (3) is converted to two-dimensional linear function, by setting up the Delanay triangulation network or utilizing redundant network to build more complicated annexation strengthening and treat the constraint of resolving system of equations, the Coherent Targets that utilizes contiguous rule that all distances are met to atmosphere correlation distance couples together, after the mutual deviation having solved between consecutive point, solve height value and the rate of deformation field of each target with respect to reference point by least square or average weighted method;

When enough hour of the interferogram time interval and Space Baseline, distortion phase place can be ignored, and two dimensional model formula (3) can be transferred to one-dimensional model formula (6), and then resolved and obtain DEM;

$\mathrm{\&Delta;}{\mathrm{\&phi;}}_{i,j}^k=C_B\&CenterDot;B^{\left(k\right)}\&CenterDot;{\mathrm{\&Delta;\&epsiv;}}_{i,j}+{\mathrm{\&alpha;}}^{\left(k\right)}+n^{\left(k\right)}---\left(6\right)$

(1) choose ultra-short Time and Space Baseline interference image to data set

According to above-mentioned theory, according to the relation of the height of buildings, interferogram baseline and landform phase place fringe density, choose the average depth of building in study area and retrain as elevation, determine that the interference image needing while calculating object height journey is primitively to sequence; Access time baseline and Space Baseline all shorter interference image to sequence as the data set initially solving, set initial time and Space Baseline threshold value and be respectively T ₁and B ₁to the interference image of data centralization to interfering processing, carry out landform phase correction without outside DEM, the method that comprehensive utilization Fast Fourier Transform (FFT) estimation and fitting of a polynomial are estimated is removed orbit error and the tendency interference fringe in interferogram, generates ultrashort space-time base line interference diagram data collection, coherence map data set and intensity map data set; Each interferogram is carried out to phase unwrapping, according to one-dimensional model repeatedly iteration correction vertical error phase estimation obtain original DEM; It should be noted that, operations all before geocoding are all carried out under radar fix system;

(2) utilize point target recognition methods to extract Coherent Targets candidate point

For the interferogram sequence of above-mentioned all generations, comprehensively adopting amplitude dispersion index is that Amplitude Dispersion Index and coefficient of coherence coherence screen Coherent Targets candidate point;

The computing formula of amplitude dispersion index is:

$D_A=\frac{{\mathrm{\&sigma;}}_A}{m_A}---\left(7\right)$

Wherein, σ _aand m _abe respectively standard deviation and the average of pixel amplitudes value; A given appropriate threshold value d _abe defined as Coherent Targets candidate point lower than the pixel of threshold value;

The Coherence Estimation formula of radar interference phase diagram is:

$\widetilde{\mathrm{\&gamma;}}=\left|\frac{\frac{1}{N}{\mathrm{\&Sigma;}}_{i=0}^N{M}_i{S}_i^{*}}{\frac{1}{N}{\mathrm{\&Sigma;}}_{i=0}^N{M}_i{M}_i^{*}\frac{1}{N}{\mathrm{\&Sigma;}}_{i=0}^N{S}_i{S}_i^{*}}\right|---\left(8\right)$

Coefficient of coherence sequence γ according to each pixel point in coherence map _iwith given coefficient of coherence threshold value if mean so this pixel is defined as to Coherent Targets candidate point;

(3) iteration correction vertical error phase estimation atural object elevation repeatedly

For each Coherent Targets point, utilize ultrashort space-time base line interference diagram data collection and according to repeatedly iteration correction vertical error phase place of one-dimensional model formula (6) Suo Shu, the vertical error correction phase place of repeatedly estimating is added, obtain the elevation at Coherent Targets point place, then, point vector is converted to raster data, generates the original DEM of study area;

Step 3: choose InSAR time series analysis interference image to data set and extract Coherent Targets candidate point

All interference images that Space Baseline and time basis are no more than to a certain setting threshold are to interfering processing, and the original elevation artificially generated terrain phase place that step 2 is generated, for the elevation phase compensation of all interferograms, generates differential interferometry graphic sequence; For the orbit error occurring in single differential interferometry figure and tendency interference fringe, the method that comprehensive utilization Fast Fourier Transform (FFT) estimation and fitting of a polynomial are estimated is removed, final differential interferometry diagram data collection and corresponding coherence map data set and the intensity map data set generating for InSAR time series analysis; To the differential interferometry graphic sequence of above-mentioned all generations, comprehensively adopt amplitude dispersion index and coefficient of coherence to screen Coherent Targets candidate point, reduce Coherent Targets quantity redundancy under high resolution SAR condition;

Step 4: elevation phase correction and linear deformation component are estimated

Solution twines each the width differential interferometry figure in InSAR time series analysis differential interferometry graphic sequence; Access time and Space Baseline threshold value are respectively T ₂and B ₂interference image to data set, re-use one-dimensional model secondary and estimate elevation residual value dh _ori; By above-mentioned one-dimensional model gained elevation residual value dh _orias the initial elevation of deformation time series analysis two-dimensional parameter estimation model, the step-length using Δ T, Δ B and Δ dh as time, Space Baseline and elevation correction, by T ₂+ Δ T and B ₂interferogram in+Δ B participates in the sequence of two-dimensional parameter estimation, sets initial maximum elevation correction threshold DH _{max_ori}, and estimate for the first time elevation modified value and rate of deformation; Progressively increase the interferogram quantity for iterative estimation with step delta T and Δ B respectively, and progressively reduce maximum elevation correction threshold with step delta dh, progressively revise elevation estimated value and rate of deformation by iterative processing repeatedly; All participate in calculating until meet all differential interferometry figure of time and Space Baseline and maximum elevation correction threshold, estimate final elevation correction and rate of deformation; Be the estimation to removing height value because successive iteration solves the elevation correction obtaining, thereby final atural object elevation estimated value is successive iteration modified value sum, that is:

$H=\underset{i=1}{\overset{u}{\mathrm{\&Sigma;}}}{\mathrm{\&epsiv;}}^i---\left(9\right)$

Atural object elevation estimated value when H is iteration termination in formula, u is total iterations;

Step 5: the deformation sequence estimation of Coherent Targets

For thering is remarkable non-linear deformation process, still need residual phase to carry out more complicated processing, to extract non-linear deformation quantity; The prerequisite of processing is still the spatio-temporal frequency characteristic based on out of phase, and the residual phase of removing from differential phase after elevation and linear velocity is:

${\mathrm{\&Delta;w}}_{i,j}^k={\mathrm{\&mu;}}_{i,j}^k+a_{i,j}^k+n_{i,j}^k---\left(10\right)$

Owing to calculating complete cycle phase place, thereby residual phase is phase place main value, and its size meets:

$\left|{\mathrm{\&Delta;w}}_{i,j}^k\right|<\mathrm{\&pi;}---\left(11\right)$

Here first to realize the filtering processing to residual phase in single differential interferometry figure; Because atmosphere shows space low frequency characteristic, and the spatial variations scope of nonlinear deformation is less, and relatively atmospheric phase shows as high pass characteristic, thereby, residual phase bitmap is carried out low pass spatial filtering and is further weakened the impact of atmosphere; The non-linear deformation quantity solving under short baseline set condition, is different from the processing policy of PSInSAR, directly atmosphere is not estimated, but is weakened the impact of atmosphere by filtering; Residual phase is after the high-pass filtering of space, and atmosphere component weakens, and now, further processing is to solve corresponding with radar data not nonlinear deformation amount in the same time, comprises noise effect; According to the relation of major-minor image, the residual phase after phase unwrapping is decomposed into:

${\mathrm{\&Delta;w}}_i^k={\mathrm{\&Delta;w}}_i^p-{\mathrm{\&Delta;w}}_i^q,\&ForAll;k=1,\&CenterDot;\&CenterDot;\&CenterDot;,N---\left(12\right)$

In formula, p and q represent to generate k and open the acquisition time of the major-minor image of differential interferometry figure; Owing to adopting short baseline principle, solving each moment when corresponding residual phase, there will be rank defect system of equations, the way addressing this problem be exactly singular value decomposition method be SVD; Thus, on the basis of svd, M is opened to residual phase and carries out time domain low-pass filtering treatment, to extract final non-linear deformation quantity, be expressed as:

${\left(d_i^p\right)}_{\mathrm{NL}\_\mathrm{est}}=\frac{\mathrm{\&lambda;}}{4\mathrm{\&pi;}}\&CenterDot;{\left\{\right[{\mathrm{\&Delta;w}}_i^p{]}_{\mathrm{HP}\_\mathrm{space}}\}}_{\mathrm{TP}\_\mathrm{time}}---\left(13\right)$

Solving after non-linear deformation quantity, the deformation sequence that each Coherent Targets is corresponding is:

$d^p=v_{\mathrm{est}}\&CenterDot;(t^p-t^1)+{\left(d_i^p\right)}_{\mathrm{NL}\_\mathrm{est}},k=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,M---\left(14\right)$

Wherein, v _estobtain Linear deformation rate for resolving.