CN109375222A - A kind of synthetic aperture radar interferometry ionosphere phase estimation and compensation method - Google Patents

A kind of synthetic aperture radar interferometry ionosphere phase estimation and compensation method Download PDF

Info

Publication number
CN109375222A
CN109375222A CN201811542725.8A CN201811542725A CN109375222A CN 109375222 A CN109375222 A CN 109375222A CN 201811542725 A CN201811542725 A CN 201811542725A CN 109375222 A CN109375222 A CN 109375222A
Authority
CN
China
Prior art keywords
phase
image
subpictures
resampling
interferometry
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201811542725.8A
Other languages
Chinese (zh)
Other versions
CN109375222B (en
Inventor
刘斌
葛大庆
李曼
张玲
王艳
郭小方
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Aero Geophysical Survey & Remote Sensing Center For Land And Resources
Original Assignee
China Aero Geophysical Survey & Remote Sensing Center For Land And Resources
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China Aero Geophysical Survey & Remote Sensing Center For Land And Resources filed Critical China Aero Geophysical Survey & Remote Sensing Center For Land And Resources
Priority to CN201811542725.8A priority Critical patent/CN109375222B/en
Publication of CN109375222A publication Critical patent/CN109375222A/en
Application granted granted Critical
Publication of CN109375222B publication Critical patent/CN109375222B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9094—Theoretical aspects
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/40—Means for monitoring or calibrating

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

A kind of synthetic aperture radar interferometry ionosphere phase estimation of the invention and compensation method, which is characterized in that specific step is as follows for this method: step 1: major-minor image accuracy registration and resampling subpictures；Step 2: look-up table optimization and again resampling subpictures；Step 3: the bandpass filtering and full bandwidth of major-minor image, subband differential interferometry map generalization；Step 4: ionosphere phase estimation；Step 5: Ground Deformation phase resolves.The present invention can obtain more accurate registration offset, and the coherence of image significantly improves, and interference line figure is more continuous.The present invention constructs front and back subband SLC image, can identify that ionospheric path postpones the influence to radar image using multiple aperture SAR interferometry method.This method constructs the upper and lower subband signal of two different center frequencies, estimates processing step by optimization ionosphere, estimates and eliminates phase error caused by the electron density variation due to ionosphere, improve interferometry precision.

Description

A kind of synthetic aperture radar interferometry ionosphere phase estimation and compensation method
Technical field
The present invention relates to synthetic aperture radar interferometry (InSAR) ionosphere phase estimation and compensation methodes, belong to conjunction At aperture SAR interferometry technical field.It can effectively weaken interfering synthetic aperture radar in the case where coherence meets Ionosphere delay phase effect in measurement process, so as to obtain more accurate Ground Deformation figure.
Background technique
When long wavelength (such as L-band) spaceborne radar (SAR) system imaging, electromagnetic wave, which passes through, postpones phase caused by ionosphere Radar imagery quality can be influenced with elevation phase aliasing.Ionosphere delay phase can be divided mainly into single order phase and second order phase, Wherein single order phase can generate group delay effect, lead to the effect of dispersion meeting that target imaging position is inaccurate, and second order phase generates Satellite-borne SAR signal pulse is caused to broaden and be distorted.For repeat track interference SAR system, major-minor image conjugate multiplication is obtained Interferometric phase may be expressed as:
Wherein, Δ φNondispThe respectively on-dispersives phase such as elevation, elevation variation, tropospheric delay, Δ φiono、Δ φMag、ΔφBendRespectively single order, second order and three rank ionosphere phases, f0It is respectively that original carrier frequency and upper and lower subband carry with f Frequently.ΔφMagIndicate the Faraday rotation effect in earth's magnetic field, Δ φBendIndicate the simple bending such as high-order earth magnetism, bending of a ray ?.Δ φ is thought by high-precision GNSS location data researchMagWith Δ φBendThe order of magnitude of this two sums is than Δ φionoIt is small It is more, for InSAR double frequency phase data, be generally only capable of for estimating Δ φNondispWith Δ φiono.Therefore, formula (1) can Simplify are as follows:
Based on this it is assumed that carrying out distance to spectrum imaging (as schemed to haplopia plural number image (SLC) using bandpass filter Shown in 1), it can be achieved that ionosphere phase estimation, to obtain more accurate Ground Deformation information.
Summary of the invention
1. purpose: the object of the present invention is to provide a kind of synthetic aperture radar interferometry ionosphere phase estimation and compensation Method, it realizes accuracy registration using azimuth segmentation spectral method iterative estimate registration offset, improves interference coherence, Traditional distance is simplified simultaneously to the process of split spectrum estimation ionosphere phase, and ionosphere delay phase can be effectively estimated.
2. technical solution: the present invention is a kind of synthetic aperture radar interferometry ionosphere phase estimation and compensation method, Specific step is as follows for this method:
Step 1: major-minor image accuracy registration and resampling subpictures
Major-minor image accuracy registration is the important step for realizing SAR interferometry, position of the same atural object on major-minor image Setting offset mainly includes that position caused by radar imaging system stability, adjacent picture elements elevation difference and surface displacement etc. is inclined It moves, for long wavelength's radar system, ionosphere generates the offset that group delay effect also results in target imaging position.In the present invention In, major-minor image accuracy registration and resampling subpictures are carried out using the look-up table and polynomial fit function of optimization.
Using distance-Doppler model (R-D model) and earth model, imageable target is established in master according to formula (3) Nonlinear Mapping relationship between image coordinate and subpictures coordinate.
Wherein, R is that sensor arrives the distance between point target, RsAnd RtThe respectively position vector of sensor and point target, λ is radar wavelength, fdFor Doppler frequency, vsAnd vtThe respectively velocity vector of sensor and point target, ReFor equatorial positions The radius of a ball, RpFor polar region radius, h is point target elevation, (xt,yt,zt) it is point target position vector under Earth central inertial reference system Component.The orientation and distance that ground point target can be calculated using formula (3) exist to imaging time to obtain point target Position on radar image.Hence, it can be determined that corresponding coordinate value of the pixel on major-minor image, can be obtained initial lookup Table.
Using the digital elevation model (DEM) under major-minor image and geographic coordinate system, is calculated according to formula (3) and initially looked into Look for table, initial look-up table is contained with main image cell coordinate at a distance under corresponding subpictures coordinate system and sat to orientation Mark.Main image is transformed under subpictures coordinate system using initial look-up table, is calculated and is turned using intensity correlation back tracking method (formula 4) Offset is registrated between the main image pixel changed under subpictures coordinate system and corresponding subpictures pixel:
Wherein, p and q respectively indicates the two width images for participating in calculating, and C expression is arranged in two width images with the i-th row j as center N row M column two windows covariance function, n and m respectively indicate row and column offset, and k and l indicate to participate in phase in window The column locations for calculating pixel are closed, when (Δ i, Δ j) take (nmax,mmax) when indicate C obtain maximum value when ranks number, can obtain Offset is registrated between main image pixel and corresponding subpictures pixel under to subpictures coordinate system.It is right using formula (5) The registration offset and subpictures cell coordinate of above-mentioned calculating carry out fitting of a polynomial:
Wherein, distance to orientation is registrated offset to Rg_off and Az_off respectively, rg and az be respectively distance to Orientation coordinate, aiAnd biParameter respectively to be estimated.Initial look-up table is optimized using the multinomial of fitting, utilization is above-mentioned The initial look-up table of optimization will be under subpictures resampling to main image geometric Framework.
Since there are errors for DEM, track condition vector and geocoding etc., only carrying out resampling with look-up table can be in weight There are a small amount of offsets between the subpictures and main image of sampling, it is therefore desirable to multinomial to deviating with orientation registration in conjunction with distance Formula carries out resampling to subpictures again.After calculating main image and the 1st resampling using intensity correlation back tracking method (formula 4) The registration offset of subpictures is registrated offset and main shadow between main image and subpictures after resampling using formula (5) As coordinate carries out fitting of a polynomial.Based on being registrated offset multinomial and optimization between main image and subpictures after resampling Under look-up table resampling subpictures to main image radar fix system.
Step 2: look-up table optimization and again resampling subpictures
Multiple aperture SAR interferometry method (MAI) carries out radar haplopia plural number image (SLC) using bandpass filter Azimuth spectrum segmentation, generates two pairs of major-minor images of forward and backward view, forward sight and backsight major-minor image is carried out interference respectively and handled Interfere line figure to forward sight and backsight, this two forward and backward view interference pattern conjugate multiplications are then generated into a width multiple aperture differential interferometry Figure, to obtain orientation phase changing capacity φMAI.L is radar antenna length, φMAIWith atural object location variation (as deformed) x Functional relation are as follows:
Influence of the ionosphere to registration accuracy can be identified using the differential interferometry phase that MAI method calculates.By orientation Split spectrum differential interferometry phase transition is orientation pixel offset deltaaz, orientation pixel offset and azimuth segmentation frequency spectrum The relationship of differential interferometry phase are as follows:
Wherein, λ is wavelength, θfAnd θbRespectively forward and backward strabismus angle.
Look-up table Optimization Steps: using the 2nd resampling subpictures in formula (6) and (7) estimation main image and step 1 it Between orientation pixel offset, the orientation pixel offset of estimation is added with the optimization look-up table in step 1, into one Then step optimization look-up table utilizes being registrated between new look-up table and the step 1 main image estimated and subpictures after resampling Deviate multinomial resampling subpictures again.The ionosphere of the subpictures of main image and resampling is reevaluated using MAI method It influences, if above-mentioned look-up table Optimization Steps are repeated there is also obvious non-zero phase in the image generated, until being not present Until apparent non-zero phase.
Step 3: the bandpass filtering and full bandwidth of major-minor image, subband differential interferometry map generalization
Major-minor SLC image M and S generates subband SLC image on respective major-minor using band-pass filter in orientation Mh、ShAnd subband SLC image M under major-minorl、Sl.The normalization centre frequency and bandwidth of upper subband bandpass filter be all 0.333；The normalization centre frequency and bandwidth of lower subband bandpass filter are respectively 0.666 and 0.333.Radar receiver passes through It is mixed with reference signal (carrier wave), the echo that microwave frequency obtains is converted into base usually at the centre frequency of radar frequency spectrum Band, bandwidth needed for this method energy minimization expresses SLC information.Bandpass filtering obtain upper and lower subband SLC image be not also Base band, this is because the deviation of new carrier frequency and original carrier frequency can cause phase modulation.Therefore, it is necessary to be solved in time-domain It adjusts, in oblique distance upwards multiplied by demodulation phase, to obtain the SLC image with new upper and lower carrier frequency.
Full bandwidth, subband interference processing pass through corresponding major-minor image conjugate multiplication φ=MS*It obtains, utilizes landform height Differential interferometry phase delta phi is obtained after journey (DEM) simulation interferometric phase compensation landform phase.Wherein, upper and lower subband differential interferometry Processing is completed using following steps:
Offset multinomial and step 2 life are registrated using between the main image and subpictures after resampling of step 1 estimation At optimization map table, respectively to upper and lower subband subpictures ShAnd SlResampling is carried out, to the upper subband major-minor image after registration Interference processing is carried out respectively with lower subband major-minor image, and landform phase is eliminated using corresponding altitude data, it can according to formula (2) Obtain upper and lower subband differential interferometry phase delta phihWith Δ φl:
Wherein, fhAnd flRespectively upper and lower subband carrier frequency.
Step 4: ionosphere phase estimation
Full bandwidth, upper and lower subband differential interferometry phase delta phi0、ΔφhWith Δ φlIt is wound around, in estimation ionosphere phase Solution is needed to twine before position.Before solution twines, the present invention filters differential interferometry figure using Goldstein-Werner filtering method Wave can reduce noise and interpolation Low coherence region.Space solution is carried out to differential interferometry figure using minimum cost flow algorithm (MCF) It twines, it is generally the case that phase signal and coherence's spatial variations are significant, can twine processing using iterative solution: at the beginning of differential interferometry figure Filtering begin using strong filtering and solution twines, obtains the space low-frequency component figure of a width differential interferometry figureTo the 1st differential interferometry Residual plotRelatively filter by force and solution is extortedIfWhole phase does not control in (- π, π), after Continue to the 2nd differential interferometry residual plotRelatively filter by force and solution is extortedIt repeats the above process, directly ExtremelyBasic control, will in (- π, π)ExtremelySuperposition is that the solution of Δ φ twines phase.Solution twines phase and continues using difference Filter window and the moving window filter of weight be iterated space filtering.By resolving equation group (8), can be ionized Layer estimation phase
Formula (9) is carried out conversion and calculates ionosphere phase by the present invention, to reduce processing step, reduces phase noise reconciliation Twine error.
For upper and lower subband carrier frequency fh=f0+ Δ f/2 and fl=f0Δ f/2, due to carrier shift amount Δ f < < f0, because This and b can simplify are as follows:
Formula (9) can be rewritten are as follows:
Step 5: Ground Deformation phase resolves
Differential interferometry phase delta phi is deformation phase delta phidefo, ionosphere phase delta phiionoAnd track landform convection current The residual phases Δ φ such as layerotherSummation.For long wavelength's radar data, Δ φotherIt is relatively small negligible, therefore According to the ionosphere phase that formula (10) are estimated, deformation phase delta phidefoEstimated valueAre as follows:
3. advantage and effect: the present invention provides a kind of synthetic aperture radar interferometry ionosphere phase estimation and compensation Method, advantage are as follows:
(1) this method deviates the look-up table accuracy registration major-minor image and resampling pair shadow of multinomial and optimization using registration Picture.Orientation gradient in ionospheric path delay can lead to the offset of orientation position, and L-band SAR image is up to several pictures Member offset is only capable of obtaining whole registration offset based on registration offset multinomial, can not compensate ionospheric path delay gradient Caused local effect, this method can obtain more accurate registration offset, and the coherence of image significantly improves, and interfere line Scheme more continuous.
(2) this method constructs front and back subband SLC image, can identify ionization using multiple aperture SAR interferometry method Influence of the layer path delay to radar image.
(3) this method constructs the upper and lower subband signal of two different center frequencies, passes through optimization ionosphere estimation processing Step is estimated and is eliminated phase error caused by the electron density variation due to ionosphere, improves interferometry precision.
Detailed description of the invention
Fig. 1 (a) carrier frequency f0And bandwidth BwFull bandwidth distance to spectrum diagram, the point target of range sensor R is corresponding Base band SLC data are as follows:Wherein T and τ is respectively pulse length and apart from point target Delay time 2R/c, c is the light velocity, sinc function is defined as: sinc (x)=sin (π x)/π x.(b) subband carrier frequency f onhAnd band It is wideBandpass filtering distance to spectrum diagram, the corresponding band logical SLC data of the point target of range sensor R are as follows:Wherein ThFor pulse length；Lower subband carrier frequency flAnd bandwidth's Bandpass filtering distance is to spectrum diagram, the corresponding band logical SLC data of the point target of range sensor R are as follows:Wherein TlFor pulse length.(c) subband carrier frequency f onhAnd bandwidth Base band distance to spectrum diagram, the corresponding base band SLC data of the point target of range sensor R are as follows:Lower subband carrier frequency flAnd bandwidthBase band distance to spectrum diagram, distance The corresponding base band SLC data of the point target of sensor R are as follows:
Fig. 2 flow diagram of the present invention.
Fig. 3 (a) 2007/11/21SLC image ionospheric effect tentatively identifies figure.
Fig. 3 (b) 2008/01/06SLC image ionospheric effect tentatively identifies figure.
Wherein, white belt-like strip indicates non-zero registration offset, and it is significant to show that the scape data ionospheric effect influences.
The orientation phase diagram that Fig. 4 (a) tradition registration resampling obtains.
Wherein, grey belt-like strip indicates that non-zero phase, ionospheric effect influence significant.
The orientation phase diagram that Fig. 4 (b) present invention registration resampling obtains.
Coherence's figure that Fig. 5 (a) tradition registration resampling obtains.
Wherein, grey belt-like strip indicates to be reduced by ionosphere effect coherence.
The coherence map that Fig. 5 (b) present invention registration resampling obtains.
The differential interferometry figure that Fig. 6 (a) present invention obtains.
The upper subband differential interferometry figure that Fig. 6 (b) present invention obtains.
The lower subband differential interferometry figure that Fig. 6 (c) present invention obtains.
The ionosphere phase diagram that Fig. 6 (d) present invention estimates.
Differential interferometry figure after the removal ionosphere phase that Fig. 7 present invention obtains.
Specific embodiment
By taking ALOS-1PALSAR monitors Djibouti port region Ground Deformation as an example, illustrate the present invention in practical engineering applications Concrete operation step.As shown in Fig. 2, the present invention is a kind of synthetic aperture radar interferometry ionosphere phase estimation and compensation Method, specific step is as follows for this method:
Step 1: the resampling of major-minor image accuracy registration and subpictures
(1) SAR data is chosen
The 2 scape ALOS-1PALSAR radar datas in covering Djibouti port region are chosen, acquisition time is November 21 in 2007 Day (2007/11/21) and January 6 (2008/01/06) in 2008, since this area is relatively stable, and time interval is shorter, can Think that the interferometric phase that differential interferometry processing measures should be in principle 0.PALSAR is the master that Japanese ALOS-1 satellite carries Dynamic formula microwave remote sensor, has 3 kinds of observation modes such as band, polarization and scanning, and revisiting period 46 days.Radar key data parameter As shown in table 1.
Table 1: ALOS-1PALSAR radar major parameter table is selected
Carrier frequency 1.2700000e+009Hz
Carrier wavelength 23.5cm
Pulse bandwidth 2.8000000e+007Hz
Incidence angle 38.8°
Distance is to Pixel size 4.684862m
Orientation Pixel size 3.181560m
Distance is to sample rate 1.2000000e+008Hz
Radar scanning mode Band pattern
Data type Haplopia plural number (SLC)
(2) ionospheric effect detects
The SAR data that 2007/11/21 is obtained generates forward and backward view SLC image using bandpass filter, utilizes formula 4 Obtain the registration offset between forward and backward view SLC image in orientation；Same step obtains 2008/01/06 forward and backward view Registration offset between SLC image in orientation.As shown in Fig. 3 (a) and Fig. 3 (b), the image of acquisition on January 6th, 2008 Registration offset in orientation shows the scape image by apparent ionosphere there are apparent ribbon non zero-offset amount Effects.
(3) major-minor image accuracy registration and resampling subpictures
Using 2007/11/21 image as main image, the Parameter File of 2007/11/21 and 2008/01/06 liang of scape image is utilized (containing orbit information) and dem data calculate initial look-up table using formula (3), obtain and 2007/11/21 image picture element coordinate phase The distance of corresponding 2008/01/06 image picture element to orientation coordinate.2007/11/21 image is turned using initial look-up table It changes under 2008/01/06 image radar fix system, is transformed under 2008/01/06 coordinate systems in image using the calculating of formula 4 Be registrated offset between 2007/11/21 image and 2008/01/06 image, using 5 pairs of formula estimate registration offset with 2008/01/06 image coordinate carries out fitting of a polynomial, is then optimized using the multinomial of fitting to initial look-up table, base In 2008/01/06 image of look-up table resampling of optimization.
The offset that is registrated of 2007/11/21 image and 2008/01/06 image after resampling, benefit are calculated using formula 4 Fitting of a polynomial is carried out to the registration offset of calculating and 2007/11/21 image coordinate with formula (5).Based on 2007/11/21 Image is adopted again again with the look-up table for being registrated offset polynomial fitting and optimization of 2008/01/06 image after resampling 2008/01/06 image of sample.(present case is estimated using 3 parametric polynomials: the parameter to be estimated of calculating is respectively a0=0.25709, a1 =2.48639e-005, a2=-3.28442e-005, b0=0.0085, b1=-6.47649e-005, b2=1.75214e- 005；Registration accuracy: distance is to 0.0118, orientation 0.23)
Step 2: look-up table optimization and subpictures resampling again
Prepare 2008/01/06 image in 2007/11/21 image and step 1 after resampling, is calculated using formula (6) There is apparent non-zero ribbon grain in the figure as shown in Fig. 4 (a) in phase changing capacity, shows to ionize in InSAR interference line figure Layer influences seriously, to need further accuracy registration.
Above-mentioned calculating phase changing capacity is converted into orientation pixel offset, and the step 1 that is added to using formula (7) In optimization initial look-up table in, advanced optimize look-up table, using the look-up table of generation and 2007/11/21 image with The registration offset polynomial fitting of 2008/01/06 image after resampling is to the resampling again of 2008/01/06 image.It will 2008/01/06 image of 2007/11/21 image and resampling estimates ionosphere effect using formula (6), such as Fig. 4 (b) institute Show, apparent non-zero phase is substantially not present, shows that accuracy registration is met the requirements.
It is coherence's figure that step 1 resampling obtains shown in Fig. 5 (a), there are apparent Low coherence region bands；Fig. 5 It (b) is coherence's figure that step two re-sampling obtains shown in, Low coherence region band disappears, and coherence significantly improves.
Step 3: the bandpass filtering and full bandwidth of major-minor image, subband differential interferometry map generalization
2008/01/06 image of 2007/11/21 image and step 2 resampling is subjected to interference processing, utilizes DEM Differential interferometry figure is obtained after simulation interferometric phase compensation landform phase, as shown in Fig. 6 (a).
The upper and lower subband SLC image of 2007/11/21 and 2008/01/06 image is generated using band-pass filter, is carried Frequency is respectively 1.2793324e+009Hz and 1.2606676e+009Hz, bandwidth 9335200.0Hz.Utilize step 1 estimation 2007/11/21 image and 2008/01/06 image after resampling be registrated offset polynomial fitting and step 2 calculate it is excellent Change look-up table, resampling is carried out to upper and lower 2008/01/06 image of subband respectively, to the upper and lower subband major-minor image after registration Interference processing is carried out respectively, and landform phase is eliminated using corresponding DEM simulation interferometric phase, obtains upper and lower subband differential interferometry Figure, as shown in Fig. 6 (b) and Fig. 6 (c).
Step 6: ionosphere phase estimation
4 filters are carried out to full bandwidth and upper and lower subband differential interferometry figure using Goldstein-Werner filtering method Wave, coherence's threshold value are respectively set to 0.25,0.35,0.45,0.5, using minimum cost flow algorithm (MCF) to filtered difference Divide interference pattern to carry out space solution to twine, solution is twined phase and is filtered using the moving window filter of 64 filter windows.From original difference Above-mentioned filtered solution is subtracted in point interference pattern and twines phase, continues to carry out residual volume 4 filtering, coherence's threshold value is respectively set It is 0.25,0.35,0.45,0.5, space is carried out to filtered remaining difference interference pattern using minimum cost flow algorithm (MCF) Solution twines, and solution is twined phase and is filtered using the moving window filters of 32 filter windows, after filtering phase control substantially (- π, In π), twined without continuing filter solution.It is that final solution twines phase that solution after this 2 times filtering, which is twined Phase Stacking,.
It is -34.019481 that a using formula (11) calculating, which is 0.499973, b, then estimates electricity using formula (10) Absciss layer phase, as shown in Fig. 6 (d).
Step 7: Ground Deformation phase resolves
Full bandwidth differential interferometry figure is filtered using Goldstein-Werner filtering method, the setting of coherence's threshold value It is 0.4, minimum cost flow algorithm (MCF) carries out space solution to filtered differential interferometry figure and twines, and subtracts the electricity of step 5 estimation Absciss layer phase is Ground Deformation phase, as shown in Figure 7.

Claims (7)

1. a kind of synthetic aperture radar interferometry ionosphere phase estimation and compensation method, which is characterized in that this method is specific Steps are as follows:
Step 1: major-minor image accuracy registration and resampling subpictures
Major-minor image accuracy registration and resampling subpictures are carried out using the look-up table and polynomial fit function of optimization；
Using distance-Doppler model, that is, R-D model and earth model, imageable target is established in main image according to formula (3) Nonlinear Mapping relationship between coordinate and subpictures coordinate；
Wherein, R is that sensor arrives the distance between point target, RsAnd RtRespectively the position vector of sensor and point target, λ are Radar wavelength, fdFor Doppler frequency, vsAnd vtThe respectively velocity vector of sensor and point target, ReFor the equatorial positions earth Radius, RpFor polar region radius, h is point target elevation, (xt,yt,zt) it is point target position vector under Earth central inertial reference system Component；The orientation and distance for calculating ground point target using formula (3) are to imaging time, to obtain point target in radar shadow As upper position；Therefore, corresponding coordinate value of the pixel on major-minor image, i.e. acquisition initial look-up table be can determine that；
Using the digital complex demodulation under major-minor image and geographic coordinate system, initial look-up table is calculated according to formula (3), Initial look-up table contain with main image cell coordinate at a distance under corresponding subpictures coordinate system to orientation coordinate；It will Main image is transformed under subpictures coordinate system using initial look-up table, is transformed into using intensity correlation back tracking method, that is, calculating of formula 4 Offset is registrated between main image pixel and corresponding subpictures pixel under subpictures coordinate system:
Wherein, p and q respectively indicates the two width images for participating in calculating, and C indicates to arrange the N for center with the i-th row j in two width images The covariance function of two windows of row M column, n and m respectively indicate row and column offset, and k and l indicate to participate in correlometer in window The column locations for calculating pixel, when (Δ i, Δ j) take (nmax,mmax) when indicate C obtain maximum value when ranks number, obtain subpictures Offset is registrated between main image pixel and corresponding subpictures pixel under coordinate system；Using formula (5) to above-mentioned calculating Registration offset and subpictures cell coordinate carry out fitting of a polynomial:
Wherein, distance to orientation is registrated offset to Rg_off and Az_off respectively, rg and az be respectively distance to and orientation To coordinate, aiAnd biParameter respectively to be estimated；Initial look-up table is optimized using the multinomial of fitting, utilizes above-mentioned optimization Initial look-up table will be under subpictures resampling to main image geometric Framework；
Step 2: look-up table optimization and again resampling subpictures
Multiple aperture SAR interferometry method MAI carries out orientation frequency to radar haplopia plural number image SLC using bandpass filter Spectrum segmentation, generates two pairs of major-minor images of forward and backward view, by forward sight and backsight major-minor image carry out respectively interference handle to obtain forward sight with Backsight interferes line figure, this two forward and backward view interference pattern conjugate multiplications is then generated a width multiple aperture differential interferometry figure, to obtain Obtain orientation phase changing capacity φMAI；L is radar antenna length, φMAIWith the functional relation of atural object location variation x are as follows:
The differential interferometry phase calculated using MAI method identifies influence of the ionosphere to registration accuracy；By azimuth segmentation frequency spectrum Differential interferometry phase transition is orientation pixel offset deltaaz, orientation pixel offset and azimuth segmentation spectral difference are interfered The relationship of phase are as follows:
Wherein, λ is wavelength, θfAnd θbRespectively forward and backward strabismus angle；
Step 3: the bandpass filtering and full bandwidth of major-minor image, subband differential interferometry map generalization
Major-minor SLC image M and S generates subband SLC image M on respective major-minor using band-pass filter in orientationh、Sh, And subband SLC image M under major-minorl、Sl；Radar receiver is by being mixed with reference signal, at the centre frequency of radar frequency spectrum The echo that microwave frequency obtains is converted into base band, bandwidth needed for this method energy minimization expression SLC information；Bandpass filtering The upper and lower subband SLC image obtained is not also base band, this is because the deviation of new carrier frequency and original carrier frequency can cause to modulate phase Position；Therefore, it is necessary to be demodulated in time-domain, in oblique distance upwards multiplied by demodulation phase, so that obtaining has new upper and lower load The SLC image of wave frequency rate；
Full bandwidth, subband interference processing pass through corresponding major-minor image conjugate multiplication φ=MS*It obtains, is simulated and interfered using DEM Differential interferometry phase delta phi is obtained after phase compensation landform phase；Wherein, upper and lower subband differential interferometry processing uses following steps It completes:
Using step 1 estimation main image and subpictures after resampling between be registrated offset multinomial and step 2 generation Optimize map table, respectively to upper and lower subband subpictures ShAnd SlResampling is carried out, to the upper subband major-minor image after registration under Subband major-minor image carries out interference processing respectively, eliminates landform phase using corresponding altitude data, obtains according to formula (2) Upper and lower subband differential interferometry phase delta phihWith Δ φl:
Wherein, fhAnd flRespectively upper and lower subband carrier frequency；
Step 4: ionosphere phase estimation
Full bandwidth, upper and lower subband differential interferometry phase delta phi0、ΔφhWith Δ φlBe wound around, estimation ionosphere phase it Before need solution to twine；Before solution twines, differential interferometry figure is filtered using Goldstein-Werner filtering method, reduces noise And interpolation Low coherence region；Space solution is carried out to differential interferometry figure using minimum cost flow algorithm MCF to twine；
Twine processing using iterative solution: differential interferometry figure initial filter is using strong filtering and solution twines, and obtains a width differential interferometry figure Space low-frequency component figureTo the 1st differential interferometry residual plotRelatively filter by force and solution is extortedIfWhole phase does not control in (- π, π), continues to the 2nd differential interferometry residual plotCarry out relatively strong filter Wave and solution is extortedIt repeats the above process, untilBasic control, will in (- π, π)ExtremelySuperposition is Δ φ Solution twine phase；Solution is twined phase and continues to be iterated space filter using the moving window filter of different filter window and weight Wave；By resolving equation group (8), ionosphere estimation phase is obtained
Formula (9) is subjected to conversion and calculates ionosphere phase, to reduce processing step, phase noise reconciliation is reduced and twines error；
For upper and lower subband carrier frequency fh=f0+ Δ f/2 and fl=f0Δ f/2, due to carrier shift amount Δ f < < f0, therefore a and B simplifies are as follows:
Formula (9) is rewritten are as follows:
Step 5: Ground Deformation phase resolves
Differential interferometry phase delta phi is deformation phase delta phidefo, ionosphere phase delta phiionoAnd it is track, landform, tropospheric Residual phase Δ φotherSummation；For long wavelength's radar data, Δ φotherIt is relatively small, it ignores, therefore according to public affairs The ionosphere phase of formula (10) estimation, deformation phase delta phidefoEstimated valueAre as follows:
2. a kind of synthetic aperture radar interferometry ionosphere phase estimation according to claim 1 and compensation method, Be characterized in that: since there are errors, only carrying out resampling with look-up table can exist between the subpictures and main image of resampling A small amount of offset, it is therefore desirable to be registrated offset multinomial again to subpictures progress resampling to orientation in conjunction with distance；Benefit Be registrated offset with what formula (4) calculated main image and the subpictures after the 1st resampling, using formula (5) to main image and Registration offset and main image coordinate after resampling between subpictures carry out fitting of a polynomial；After main image and resampling Under registration offset multinomial between subpictures and the look-up table resampling subpictures to main image radar fix system of optimization.
3. a kind of synthetic aperture radar interferometry ionosphere phase estimation according to claim 1 and compensation method, It is characterized in that: look-up table Optimization Steps: utilizing the 2nd resampling subpictures in formula (6) and (7) estimation main image and step 1 Between orientation pixel offset, the orientation pixel offset of estimation is added with the optimization look-up table in step 1, into Then one-step optimization look-up table utilizes matching between subpictures after new look-up table and the main image and resampling of step 1 estimation Standard deviates multinomial resampling subpictures again；The ionization of the subpictures of main image and resampling is reevaluated using MAI method Layer influences, if repeating look-up table Optimization Steps there is also apparent non-zero phase in the image generated, until there is no obvious Non-zero phase until.
4. a kind of synthetic aperture radar interferometry ionosphere phase estimation according to claim 1 and compensation method, Be characterized in that: in step 3, the normalization centre frequency and bandwidth of upper subband bandpass filter are all 0.333；Lower subband band logical The normalization centre frequency and bandwidth of filter are respectively 0.666 and 0.333.
5. a kind of synthetic aperture radar interferometry ionosphere phase estimation according to claim 1 and compensation method, It is characterized in that: carrier frequency f0And bandwidth BwFull bandwidth distance to frequency spectrum, the corresponding base band SLC data of the point target of range sensor R Are as follows:Wherein, T and τ is respectively pulse length and the delay apart from point target Time 2R/c, c are the light velocity, sinc function is defined as: sinc (x)=sin (π x)/π x.
6. a kind of synthetic aperture radar interferometry ionosphere phase estimation according to claim 1 and compensation method, It is characterized in that: upper subband carrier frequency fhAnd bandwidthBandpass filtering distance to frequency spectrum, the corresponding band of the point target of range sensor R Logical SLC data are as follows:Wherein ThFor pulse length；Lower son Band carrier frequency flAnd bandwidthBandpass filtering distance to frequency spectrum, the corresponding band logical SLC data of the point target of range sensor R are as follows:Wherein, TlFor pulse length.
7. a kind of synthetic aperture radar interferometry ionosphere phase estimation according to claim 1 and compensation method, It is characterized in that: upper subband carrier frequency fhAnd bandwidthBase band distance to frequency spectrum, the corresponding base band of the point target of range sensor R SLC data are as follows:Lower subband carrier frequency flAnd bandwidthBase band distance to frequency Spectrum, the corresponding base band SLC data of the point target of range sensor R are as follows:
CN201811542725.8A 2018-12-17 2018-12-17 Synthetic aperture radar interferometric ionosphere phase estimation and compensation method Expired - Fee Related CN109375222B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811542725.8A CN109375222B (en) 2018-12-17 2018-12-17 Synthetic aperture radar interferometric ionosphere phase estimation and compensation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811542725.8A CN109375222B (en) 2018-12-17 2018-12-17 Synthetic aperture radar interferometric ionosphere phase estimation and compensation method

Publications (2)

Publication Number Publication Date
CN109375222A true CN109375222A (en) 2019-02-22
CN109375222B CN109375222B (en) 2019-12-06

Family

ID=65374264

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811542725.8A Expired - Fee Related CN109375222B (en) 2018-12-17 2018-12-17 Synthetic aperture radar interferometric ionosphere phase estimation and compensation method

Country Status (1)

Country Link
CN (1) CN109375222B (en)

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110109102A (en) * 2019-04-04 2019-08-09 电子科技大学 A kind of method of SAR moving object detection and velocity estimation
CN111239735A (en) * 2020-02-28 2020-06-05 西南交通大学 Azimuth deformation field ionosphere correction method based on low-frequency SAR image
CN111239736A (en) * 2020-03-19 2020-06-05 中南大学 Single-baseline-based surface elevation correction method, device, equipment and storage medium
CN111366925A (en) * 2020-03-27 2020-07-03 长安大学 SAR offset two-dimensional deformation time sequence calculation method and system
CN112711021A (en) * 2020-12-08 2021-04-27 中国自然资源航空物探遥感中心 Multi-resolution InSAR (interferometric synthetic Aperture Radar) interactive interference time sequence analysis method
CN112946645A (en) * 2021-01-29 2021-06-11 北京理工大学 Unmanned aerial vehicle-mounted ultra-wideband SAR self-focusing method
CN113189586A (en) * 2021-04-02 2021-07-30 中国地质大学（武汉） Polarization phase optimization method based on PS target time-space coherent matrix
CN113406628A (en) * 2021-05-12 2021-09-17 中国科学院国家空间科学中心 Data resampling method for interference imaging altimeter
CN113624122A (en) * 2021-08-10 2021-11-09 中咨数据有限公司 Bridge deformation monitoring method fusing GNSS data and InSAR technology
CN114722344A (en) * 2022-03-11 2022-07-08 中国电子科技集团公司第二十九研究所 High-precision data splicing method based on phase compensation
CN115372964A (en) * 2022-10-26 2022-11-22 中国电子科技集团公司第十四研究所 Double-frequency multi-scale earth surface deformation measurement test system

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103645476A (en) * 2013-12-18 2014-03-19 中国国土资源航空物探遥感中心 Space-time homogeneous filter method for synthetic aperture radar differential interference figure sequences
CN103728604A (en) * 2013-11-19 2014-04-16 中国国土资源航空物探遥感中心 Broadband synthetic aperture radar sub-band interferometric data processing method
CN108303735A (en) * 2018-01-30 2018-07-20 单新建 The earthquake deformation acquisition methods of subband interferometry based on optimized parameter setting

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103728604A (en) * 2013-11-19 2014-04-16 中国国土资源航空物探遥感中心 Broadband synthetic aperture radar sub-band interferometric data processing method
CN103645476A (en) * 2013-12-18 2014-03-19 中国国土资源航空物探遥感中心 Space-time homogeneous filter method for synthetic aperture radar differential interference figure sequences
CN108303735A (en) * 2018-01-30 2018-07-20 单新建 The earthquake deformation acquisition methods of subband interferometry based on optimized parameter setting

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
BOCHEN ZHANG ET AL.: ""Mitigating Ionospheric Artifacts in Coseismic Interferogram Based on Offset Field Derived From ALOS-PALSAR Data"", 《IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING》 *
ZHEN LIU ET AL.: ""Joint Correction of Ionosphere Noise and Orbital Error in L-Band SAR Interferometry of Interseismic Deformation in Southern California"", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》 *
张哲远: ""GEO及GEO-LEO SAR成像及干涉处理研究"", 《中国优秀硕士学位论文全文数据库（电子期刊） 信息科技辑》 *

Cited By (21)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110109102A (en) * 2019-04-04 2019-08-09 电子科技大学 A kind of method of SAR moving object detection and velocity estimation
CN110109102B (en) * 2019-04-04 2022-05-03 电子科技大学 SAR moving target detection and speed estimation method
CN111239735A (en) * 2020-02-28 2020-06-05 西南交通大学 Azimuth deformation field ionosphere correction method based on low-frequency SAR image
CN111239735B (en) * 2020-02-28 2022-08-30 西南交通大学 Azimuth deformation field ionosphere correction method based on low-frequency SAR image
CN111239736B (en) * 2020-03-19 2022-02-11 中南大学 Single-baseline-based surface elevation correction method, device, equipment and storage medium
CN111239736A (en) * 2020-03-19 2020-06-05 中南大学 Single-baseline-based surface elevation correction method, device, equipment and storage medium
CN111366925A (en) * 2020-03-27 2020-07-03 长安大学 SAR offset two-dimensional deformation time sequence calculation method and system
CN111366925B (en) * 2020-03-27 2022-11-22 长安大学 SAR offset two-dimensional deformation time sequence calculation method and system
CN112711021A (en) * 2020-12-08 2021-04-27 中国自然资源航空物探遥感中心 Multi-resolution InSAR (interferometric synthetic Aperture Radar) interactive interference time sequence analysis method
CN112711021B (en) * 2020-12-08 2021-10-22 中国自然资源航空物探遥感中心 Multi-resolution InSAR (interferometric synthetic Aperture Radar) interactive interference time sequence analysis method
CN112946645A (en) * 2021-01-29 2021-06-11 北京理工大学 Unmanned aerial vehicle-mounted ultra-wideband SAR self-focusing method
CN112946645B (en) * 2021-01-29 2022-10-14 北京理工大学 Unmanned aerial vehicle-mounted ultra-wideband SAR self-focusing method
CN113189586B (en) * 2021-04-02 2022-06-21 中国地质大学（武汉） Polarization phase optimization method based on PS target time-space coherent matrix
CN113189586A (en) * 2021-04-02 2021-07-30 中国地质大学（武汉） Polarization phase optimization method based on PS target time-space coherent matrix
CN113406628A (en) * 2021-05-12 2021-09-17 中国科学院国家空间科学中心 Data resampling method for interference imaging altimeter
CN113624122A (en) * 2021-08-10 2021-11-09 中咨数据有限公司 Bridge deformation monitoring method fusing GNSS data and InSAR technology
CN113624122B (en) * 2021-08-10 2022-09-20 中咨数据有限公司 Bridge deformation monitoring method fusing GNSS data and InSAR technology
CN114722344A (en) * 2022-03-11 2022-07-08 中国电子科技集团公司第二十九研究所 High-precision data splicing method based on phase compensation
CN114722344B (en) * 2022-03-11 2023-10-20 中国电子科技集团公司第二十九研究所 High-precision data splicing method based on phase compensation
CN115372964A (en) * 2022-10-26 2022-11-22 中国电子科技集团公司第十四研究所 Double-frequency multi-scale earth surface deformation measurement test system
CN115372964B (en) * 2022-10-26 2022-12-27 中国电子科技集团公司第十四研究所 Double-frequency multi-scale earth surface deformation measurement test system

Also Published As

Publication number Publication date
CN109375222B (en) 2019-12-06

Similar Documents

Publication Publication Date Title
CN109375222A (en) A kind of synthetic aperture radar interferometry ionosphere phase estimation and compensation method
Fattahi et al. InSAR time-series estimation of the ionospheric phase delay: An extension of the split range-spectrum technique
Gomba et al. Toward operational compensation of ionospheric effects in SAR interferograms: The split-spectrum method
Fornaro et al. Three-dimensional multipass SAR focusing: Experiments with long-term spaceborne data
Kim et al. Correcting distortion of polarimetric SAR data induced by ionospheric scintillation
US4924229A (en) Phase correction system for automatic focusing of synthetic aperture radar
Hanssen Radar interferometry: data interpretation and error analysis
Liang et al. Measuring azimuth deformation with L-band ALOS-2 ScanSAR interferometry
Brcic et al. Estimation and compensation of ionospheric delay for SAR interferometry
US6919839B1 (en) Synthetic aperture radar (SAR) compensating for ionospheric distortion based upon measurement of the group delay, and associated methods
US6914553B1 (en) Synthetic aperture radar (SAR) compensating for ionospheric distortion based upon measurement of the Faraday rotation, and associated methods
US5923278A (en) Global phase unwrapping of interferograms
Stevens et al. Options for airborne interferometric SAR motion compensation
CN112711021B (en) Multi-resolution InSAR (interferometric synthetic Aperture Radar) interactive interference time sequence analysis method
Mouginot et al. Correction of the ionospheric distortion on the MARSIS surface sounding echoes
CN109116321B (en) A kind of phase filtering method and height measurement method of spaceborne interference imaging altimeter
Lee et al. Multibaseline TanDEM-X mangrove height estimation: The selection of the vertical wavenumber
Gomba et al. Bayesian data combination for the estimation of ionospheric effects in SAR interferograms
CN115032638A (en) Bistatic SAR (synthetic aperture radar) phase synchronization precision improving method based on compressed sensing
Fornaro et al. Minimum mean square error space-varying filtering of interferometric SAR data
Bamler et al. Split band interferometry versus absolute ranging with wideband SAR systems
Kim et al. SAR observation of ionosphere using range/azimuth sub-bands
Cotton Special problems in imaging
Arora et al. Ionospheric Modelling using GPS to Calibrate the MWA. II: Regional ionospheric modelling using GPS and GLONASS to estimate ionospheric gradients
CN111239735B (en) Azimuth deformation field ionosphere correction method based on low-frequency SAR image

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB02 Change of applicant information

Address after: 100083 Beijing Haidian District North Fourth Ring Road 267 Olympic Building

Applicant after: CHINA AERO GEOPHYSICAL SURVEY AND REMOTE SENSING CENTER FOR NATURAL RESOURCES

Address before: 100083 Beijing Haidian District North Fourth Ring Road 267 Olympic Building

Applicant before: China Aero Geophysical Survey & Remote Sensing Center for Land and Resources

CB02 Change of applicant information
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20191206

Termination date: 20211217

CF01 Termination of patent right due to non-payment of annual fee