CN109143234A - The InSAR filtering method and system of interferometric phase gradient are estimated based on FFT high-precision - Google Patents
The InSAR filtering method and system of interferometric phase gradient are estimated based on FFT high-precision Download PDFInfo
- Publication number
- CN109143234A CN109143234A CN201810779761.XA CN201810779761A CN109143234A CN 109143234 A CN109143234 A CN 109143234A CN 201810779761 A CN201810779761 A CN 201810779761A CN 109143234 A CN109143234 A CN 109143234A
- Authority
- CN
- China
- Prior art keywords
- interferometric phase
- phase gradient
- estimation
- interferometric
- window
- 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
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
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
The invention discloses a kind of InSAR filtering methods and system that interferometric phase gradient is estimated based on FFT high-precision, wherein this method includes carrying out 3 × 3 mean filters to original interference phase diagram first；In the estimation of interferometric phase gradient initial value, the quality coefficient of interferometric phase gradient is estimated simultaneously, quality coefficient is higher than the interferometric phase gradient initial value of threshold value as final estimated value, quality is removed lower than the initial value of threshold value, and the final interferometric phase gradient estimated value of the position is weighted and averaged to obtain using the higher interferometric phase gradient initial value of ambient quality.Finally, using the interferometric phase gradient in the interferometric phase gradient compensating interferometer phase filtering window for estimating to obtain, then mean filter is carried out, complete the whole process of InSAR phase filtering.The present invention solve faced in InSAR phase filtering filtering view number and terrain information keep between contradiction, overcome coherence it is poor when interferometric phase gradient estimation robustness difference difficult point.
Description
Technical field
The invention belongs to interference synthetic aperture radar topographic survey fields, more particularly to a kind of FFT high-precision that is based on to estimate
The InSAR filtering method and system of interferometric phase gradient.
Background technique
To the SAR complex pattern after accuracy registration to interference processing (i.e. conjugate multiplication) is carried out, multiple interferometric phase image is obtained.By
The influence of various decoherence factors, phase diagram directly affect at subsequent intervention phase unwrapping there are random interferometric phase noise
The complexity of reason and topographic precision, the target of InSAR phase filtering are reduced in major-minor SAR image interferometric phase
Random noise component, while reducing filter caused landform phase information loss as far as possible.
InSAR phase filtering algorithm is conducted extensive research both at home and abroad.Baran and Stewart et al. will
Smooth operation parameter in Goldstein filtering algorithm is revised as related to coherence factor, improves the filter of Goldstein algorithm
Wave performance.Suo Zhi is brave et al. to carry out piecemeal phase unwrapping to pretreated interferometric phase image, and compensation solution carries out again after twining phase
Multiple mean filter, then the iteration above process, the filter result precision of acquisition are higher.E.Trouv é and J.M.Nicolas et al.
A kind of multiple more view phase filtering algorithms based on the compensation of landform phase gradient are proposed, to the compensated last phase of landform phase gradient
The unanimity of samples that position filtering is effectively guaranteed in filter window, but the robustness that the algorithm is estimated in Low coherence regional slope
It is poor.J.S.Lee and M.P.Stewart et al. are based on sample statistics characteristic and the adaptively selected filtering sample of one group of directionality window
This, filtering accuracy is higher.The it is proposeds such as Wang Qingsong utilize Two-dimensional FFT and half-power THRESHOLD ESTIMATION landform phase, then dry to remnants
It relates to phase and carries out simplified Lee filtering, it is high which estimates landform phase efficiency, but deblocking size is difficult to determine, piecemeal
When estimating landform phase, the splicing trace between block and block is difficult to eliminate.Li Zhenfang et al. passes through neighborhood pixels around joint
There are the interferometric phases in the case of big registration error for coherence messages estimation.Deledalle et al. proposes same using the method for iteration
When estimation scene scatters amplitude, interferometric phase and coherence factor, but this method is also easy to produce number when major-minor SAR image amplitude is close
Value calculates mistake.
For the processing of InSAR phase filtering, in order to inhibit the random noise in interferometric phase, need to select and picture to be filtered
Element meets independent identically distributed sample, then carries out ensemble average to these samples.However it is difficult in practice by repeatedly observing
A large amount of independent same distribution sample is obtained, therefore usually replaces ensemble average using space average, is i.e. hypothesis interferometric phase image exists
Part meets spatial stationarity property and pixel to be filtered and surrounding pixel meet independent same distribution.In practice, due to landform or atural object
In regional area, there may be acute variation, pixel to be filtered and surrounding pixels not ensure that and meet independent same distribution, therefore
The filtering of InSAR interferometric phase faces the contradiction between sufficient filtered samples and terrain information holding.Current interferometric phase ladder
Degree estimation method perhaps lacks robustness or operand is big and time-consuming, and practicability is poor.
Summary of the invention
Technical problem solved by the present invention is having overcome the deficiencies of the prior art and provide a kind of based on FFT high-precision estimation
The InSAR filtering method and system of interferometric phase gradient, solve the filtering view number faced in InSAR phase filtering and landform is believed
Breath keep between contradiction, overcome coherence it is poor when the interferometric phase gradient estimation difficult points such as robustness difference, for high-precision
The InSAR measurement of higher degree provides support, has important application value.
The object of the invention is achieved by the following technical programs: according to an aspect of the invention, there is provided a kind of base
The InSAR filtering method of interferometric phase gradient is estimated in FFT high-precision, and described method includes following steps: step (1): to matching
Major-minor SAR complex pattern after standard carries out interference processing, generates initially interferometric phase image s again1(x, y), then to s1(x, y) progress 3 ×
3 mean filters obtain coarse filtration wave and answer interferometric phase image s2(x,y)；Wherein, x and y is respectively initially to answer in interferometric phase image arbitrarily
The two-dimensional coordinate of pixel；Step (2): interferometric phase image s is answered in coarse filtration wave2Setting interferometric phase gradient estimates window in (x, y),
Estimate to obtain initial interference phase gradient estimation figure based on FFTWithEstimate Quality MapIts
In, p and q are the two-dimensional coordinate that initial interference phase gradient estimates any pixel in figure,It is initially done for X direction
Phase gradient estimation figure is related to,Estimate to scheme for y direction initial interference phase gradient；Step (3): according to estimation matter
SpirogramReject initial interference phase gradient estimation figure In quality it is dry lower than 0.4
Phase gradient estimated value is related to, final interferometric phase gradient estimation figure is obtainedWithStep (4): using most
Whole interferometric phase gradient estimation figureWithInterferometric phase gradient in compensation filter window, to initial multiple
Interferometric phase image s1(x, y) carries out mean filter, multiple interferometric phase image s after being filtered3(x,y)。
It is right in the step (1) in the above-mentioned InSAR filtering method based on FFT high-precision estimation interferometric phase gradient
Initially answer interferometric phase image s1The method that (x, y) carries out 3 × 3 mean filters are as follows: filtering is set in interferometric phase image initially answering
Window, window size are 3 × 3, and it is average then to carry out plural number to pixel in window.
In the above-mentioned InSAR filtering method based on FFT high-precision estimation interferometric phase gradient, coarse filtration wave answers interferometric phase image
s2(x, y) isN is sum of all pixels in filter window, -1≤lx,ly≤1
For two integer variables, with lxAnd lyValue is different, (x+lx,y+ly) indicate to take (x, y) as the difference in the filter window of center
Position, x and y are respectively the two-dimensional coordinate of any pixel in initial interferometric phase image again.
In the above-mentioned InSAR filtering method based on FFT high-precision estimation interferometric phase gradient, the step (2) is specifically wrapped
Include following steps:
(2.1) interferometric phase image s is answered in coarse filtration wave2Setting interferometric phase gradient estimates that window, window size are in (x, y)
Nf×Nf, window center is located at (x, y), and the coarse filtration wave in window interferes diagram data to be g (m, n) again；
(2.2) zero padding is first carried out to g (m, n) to operate to obtain gz(m, n), then to gz(m, n) carries out two-dimentional fast Fourier
Transformation obtains 2-d spectrum Gz(u,v)；
(2.3) 2-d spectrum GzProminent frequency coordinate is (u in (u, v)max,vmax), calculate fx、fyAnd C:
fx=2 π (umax-Nz/2)/Nz
fy=2 π (vmax-Nz/2)/Nz
Wherein, | | it is operated for calculated complex modulus value, NzThe g operated for zero paddingzThe two-dimensional of (m, n), i.e. gz
The data scale of (m, n) is Nz×Nz, Gz(umax,vmax) it is to gz(m, n) carries out two-dimensional fast fourier transform and obtains two-dimentional frequency
Compose GzThe frequency spectrum of power maximum, u in (u, v)maxAnd vmaxRespectively 2-d spectrum GzThe two dimension frequency of power maximum in (u, v)
Rate value, fxAnd fyThe X direction and y direction initial interference phase of all pixels respectively in interferometric phase gradient estimation window
Potential gradient estimated value, C is that interferometric phase gradient estimates that the interferometric phase gradient of all pixels in window estimates quality coefficient, by fx
And fyIt is respectively stored in initial interference phase gradient estimation figureWithInPosition
It sets, C is stored in estimation Quality MapInPosition,For the operation that rounds up；
(2.4) position for changing interferometric phase gradient estimation window, so that new window center is located at (x+Nf, y) or (x,
y+Nf), it repeats step (2.1) and arrives (2.3), until completing the interferometric phase that estimation coarse filtration wave answers each pixel in interferometric phase image
Gradient and quality coefficient obtain initial interference phase gradient estimation figureWithAnd estimation Quality Map
In the above-mentioned InSAR filtering method based on FFT high-precision estimation interferometric phase gradient, the step (3) is specifically wrapped
Include following steps:
(3.1) it will estimate Quality MapMiddle all elements are normalized；
(3.2) quality threshold is setFor normalizing Quality MapIn be higher than quality threshold member
Element then determines initial interference phase gradient estimation figureWithEffectively, final interferometric phase gradient estimation figure
(3.3) for normalizing Quality MapIn each be lower than quality threshold element, then determine initial interference
Phase gradient estimation figureWithIn vain, interferometric phase gradient is reevaluated, final interferometric phase
Gradient estimation figure are as follows:
Wherein,The weighting weight of corresponding element in figure, k are estimated for initial interference phase gradientpAnd kqFor integer, table
Two-dimensional coordinate when showing the estimation of final interferometric phase gradient relative to (p, q),For in estimation Quality Map
(p+kp,q+kq) at estimation quality coefficient.
In the above-mentioned InSAR filtering method based on FFT high-precision estimation interferometric phase gradient, in the step (4), filter
Multiple interferometric phase image s after wave3(x, y) are as follows:
Wherein, it is sum of all pixels in filter window that window size, which is 7 × 7, N, and x and y are respectively initially to answer interferometric phase image
In any pixel two-dimensional coordinate, f'xInitially to answer pixel s in interferometric phase image1(x, y) corresponding X direction is finally interfered
Phase gradient estimated value, f'yInitially to answer pixel s in interferometric phase image1The final interferometric phase ladder of (x, y) corresponding y direction
Spend estimated value.
In the above-mentioned InSAR filtering method based on FFT high-precision estimation interferometric phase gradient, in step (2.1), window
In interfere the formula of diagram data g (m, n) as follows again:
G (m, n)=s2(x+m-Nf/2,y+n-Nf/ 2),
Wherein, 1≤m, n≤NfFor positive integer, m and n indicate that interferometric phase gradient estimates the two-dimensional coordinate of pixel in window.
In the above-mentioned InSAR filtering method based on FFT high-precision estimation interferometric phase gradient, in step (2.2), pass through
G after zero paddingz(m, n) points are Nz×Nz, Nz=128, gz(m, n) is indicated are as follows:
In the above-mentioned InSAR filtering method based on FFT high-precision estimation interferometric phase gradient, in step (3.1), it will estimate
Count Quality MapMiddle all elements, which are normalized, to be obtained by the following formula:
Wherein, max () expression is maximized operation.
According to another aspect of the present invention, a kind of InSAR that interferometric phase gradient is estimated based on FFT high-precision is additionally provided
Filtering system, comprising: the first module generates initial multiple interference for carrying out interference processing to the major-minor SAR complex pattern after registration
Phase diagram s1(x, y), then to s1(x, y) progress mean filter obtains coarse filtration wave and answers interferometric phase image s2(x,y)；Wherein, x and y
For the two-dimensional coordinate for initially answering any pixel in interferometric phase image；Second module, for answering interferometric phase image s in coarse filtration wave2(x,
Y) setting interferometric phase gradient estimates window in, estimates to obtain initial interference phase gradient estimation figure (two dimension interference phase based on FFT
Potential gradient initial value)WithEstimate Quality Map (quality coefficient)Wherein, p and q is initial
Interferometric phase gradient estimates the two-dimensional coordinate of any pixel in figure,For the estimation of X direction initial interference phase gradient
Figure,Estimate to scheme for y direction initial interference phase gradient；Third module, for according to estimation Quality Map (quality
Coefficient)Reject initial interference phase gradient estimation figureIn the lower interference of quality
Phase gradient estimated value obtains final interferometric phase gradient estimation figureWith4th module, for utilizing
Final interferometric phase gradient estimation figureWithInterferometric phase gradient in compensation filter window, to initial
Multiple interferometric phase image s1(x, y) carries out mean filter, multiple interferometric phase image s after being filtered3(x,y)。
Compared with prior art, the present invention has the following advantages:
(1) the method for the present invention estimates interferometric phase gradient first, compensates local landform phase first before interferometric phase filtering
Potential gradient improves the consistency of filtered samples, solves the filtering view number faced in InSAR phase filtering and terrain information is protected
Contradiction between holding；
(2) before the estimation of interferometric phase gradient, using 3 × 3 multiple mean filter, estimation is substantially improved in method of the invention
Robustness；
(3) method of the invention is based on FFT estimation regional area two dimension interferometric phase gradient, rather than directly utilizes FFT
Phase filtering is carried out, the piecemeal stitching error of Goldstein algorithm is eliminated, estimated efficiency and robustness are high；With E.Trouv é
Et al. propose the estimation method based on maximum likelihood compare, this method based on FFT estimation part two-dimentional interferometric phase gradient,
And each window only carries out 1 estimation, rather than in such a way that tradition draws window, therefore treatment effeciency greatly improves.With
The phase filtering method that directly carried out using FFT that Goldstein et al. is proposed is compared, and this method only utilizes FFT estimation single-frequency
Interferometric phase gradient, there is no need to carry out estimation thresholding setting, robustness is greatly improved, and is additionally based on the compensation of interferometric phase gradient
Mean filter without carry out piecemeal concatenation, avoid piecemeal stitching error；
(4) method of the invention establishes estimation Quality Map, rejects low-quality initial interference phase ladder according to quality coefficient
Estimated value is spent, the filtering accuracy in Low coherence region is improved.
Detailed description of the invention
By reading the following detailed description of the preferred embodiment, various other advantages and benefits are common for this field
Technical staff will become clear.The drawings are only for the purpose of illustrating a preferred embodiment, and is not considered as to the present invention
Limitation.And throughout the drawings, the same reference numbers will be used to refer to the same parts.In the accompanying drawings:
Fig. 1 is that the present invention is based on the InSAR filtering method flow charts that FFT high-precision estimates interferometric phase gradient；
Fig. 2 is N × N mean filter schematic diagram of the invention (N is odd number)；
Fig. 3 is the schematic diagram of interferometric phase gradient initial estimation of the invention, and initial and final interferometric phase gradient
The schematic diagram of estimation figure and the corresponding relationship of multiple interferometric phase image；
Fig. 4 is the main SAR data magnitude image in the Italian volcano Etna region of the invention；
Fig. 5 is the initial interferometric phase again that the major-minor SAR image of the invention to accuracy registration carries out that interference is handled
Scheme s1(x,y)；
Fig. 6 is of the invention to initial interferometric phase image s again1(x, y) carries out 3 × 3 mean filters, and obtained coarse filtration wave is multiple
Interferometric phase image s2(x,y)；
Fig. 7 answers interferometric phase image s in coarse filtration wave for the present invention2Setting interferometric phase gradient estimates window in (x, y), is based on
The X direction initial interference phase gradient estimation figure that FFT estimates；
Fig. 8 answers interferometric phase image s in coarse filtration wave for the present invention2Setting interferometric phase gradient estimates window in (x, y), is based on
The y direction initial interference phase gradient estimation figure that FFT estimates；
Fig. 9 answers interferometric phase image s in coarse filtration wave for the present invention2Setting interferometric phase gradient estimates window in (x, y), is based on
The estimation Quality Map for the initial interference phase gradient estimation figure that FFT estimates；
Figure 10 is the present invention according to estimation quality coefficient, rejects quality in X direction initial interference phase gradient estimation figure
Lower interferometric phase gradient estimated value, the final interferometric phase gradient estimation figure of obtained X direction；
Figure 11 is the present invention according to estimation quality coefficient, rejects quality in y direction initial interference phase gradient estimation figure
Lower interferometric phase gradient estimated value, the final interferometric phase gradient estimation figure of obtained y direction；
Figure 12 is the final interferometric phase gradient estimation figure that the present invention is obtained using estimation, the interference in compensation filter window
Phase gradient, to initial interferometric phase image s again1(x, y) carries out 7 × 7 mean filters, multiple interferometric phase image after being filtered.
Specific embodiment
Exemplary embodiments of the present disclosure are described in more detail below with reference to accompanying drawings.Although showing the disclosure in attached drawing
Exemplary embodiment, it being understood, however, that may be realized in various forms the disclosure without should be by embodiments set forth here
It is limited.On the contrary, these embodiments are provided to facilitate a more thoroughly understanding of the present invention, and can be by the scope of the present disclosure
It is fully disclosed to those skilled in the art.It should be noted that in the absence of conflict, embodiment in the present invention and
Feature in embodiment can be combined with each other.The present invention will be described in detail below with reference to the accompanying drawings and embodiments.
Fig. 1 is that the present invention is based on the InSAR filtering method flow charts that FFT high-precision estimates interferometric phase gradient.Such as Fig. 1 institute
Show, this method comprises the following steps:
The first step carries out interference processing to the major-minor SAR complex pattern after accuracy registration, generates initially interferometric phase image s again1
(x, y), then to s13 × 3 mean filters of (x, y) progress obtain coarse filtration wave and answer interferometric phase image s2(x, y), method are: initial
Multiple interferometric phase image is middle setting filter window having a size of 3 × 3, and window center coordinate is (x, y), and coarse filtration wave answers interferometric phase image
ForN is sum of all pixels in filter window, and x and y are respectively initial multiple dry
Relate to the two-dimensional coordinate of any pixel in phase diagram, -1≤lx,ly≤ 1 is two integer variables, with lxAnd lyValue is different, (x+
lx,y+ly) indicate to take (x, y) as the different location in the filter window of center.
Wherein, as shown in Figure 2 (N is odd number), gray circles represent in multiple interferometric phase image N × N mean filter in figure
Pixel, the process of N × N mean filter are as follows: filter window (such as 3 × 3) are set first centered on pixel " 1 ", in window
Pixel carries out plural average computation, and obtained calculated result is as filtered multiple interferometric phase measured value；Filter window is sliding
" 2 ", " 3 " ... are moved, then filter window is moved to " 9 ", with same by the filtering processing until completing the row all pixels
Process obtains the filter result of the second row data, is finally completed the filtering processing of whole picture interferometric phase image.It should be noted that figure
As the pixel of outermost edges can not carry out 3 × 3 filtering processings, can will before filtering multiple interferometric phase measured value as being tied after filtering
Fruit.
Second step answers interferometric phase image s in coarse filtration wave2Setting interferometric phase gradient estimates window in (x, y), is based on FFT
Estimate two-dimentional initial interference phase gradient estimation figureAnd estimation Quality MapMethod
It is:
2.1, it is answered in coarse filtration wave and interferometric phase gradient estimation window, window size N is set in interferometric phase imagef×Nf, pass through
Compromise considers estimated accuracy and estimated efficiency, Nf=8, interfere diagram data to be represented by again in window
G (m, n)=s2(x+m-Nf/2,y+n-Nf/2),1≤m,n≤Nf；
2.2, zero padding is first carried out to g (m, n) to operate to obtain gz(m, n), the g after zero paddingz(m, n) points are Nz×Nz, Nz
=128, gz(m, n) is represented by
Then to gz(m, n) carries out two-dimensional fast fourier transform (FFT2) and obtains 2-d spectrum Gz(u,v)；
2.3, prominent frequency coordinate is (u in 2-d spectrummax,vmax), interferometric phase gradient estimates own in window
The initial interference phase gradient estimated value table of pixel is shown as fx=2 π (umax-Nz/2)/Nz, fy=2 π (vmax-Nz/2)/Nz, interferometric phase ladder
The interferometric phase gradient estimation quality coefficient of all pixels is defined as in degree estimation window
| | it is operated for calculated complex modulus value, estimates the physical meaning of quality coefficient here are as follows: power in interferometric phase image 2-d spectrum
The ratio of the mean power of maximum value and other frequencies, the value range of quality coefficient are (0,1), the bigger expression of the estimated value of C
The estimation of interferometric phase gradient is more steady, by fxAnd fyIt is respectively stored in initial interference phase gradient estimation figureWithInC is stored in estimation Quality Map by positionIn
Position,For the operation that rounds up；
2.4, the position for changing interferometric phase gradient estimation window, so that new window center is located at (x+Nf, y) or
(x,y+Nf), step 2.1 to 2.3 is repeated, until completing the interferometric phase that estimation coarse filtration wave answers each pixel in interferometric phase image
Gradient and quality coefficient obtain initial interference phase gradient estimation figureWithWith estimation Quality MapP and q is the two-dimensional coordinate that initial interference phase gradient estimates any pixel in figure.
Fig. 3 is that interferometric phase gradient initial value estimates schematic diagram and interferometric phase gradient estimation figure and multiple interferometric phase
Corresponding relationship between figure.Left side is multiple interferometric phase image, and solid line and dashed rectangle respectively represent adjacent estimation window, right side
To estimate obtained interferometric phase gradient map or quality coefficient figure.When carrying out the estimation of interferometric phase gradient initial value, window ruler
Very little Nf×Nf, Nf=8.
Third step, according to estimation Quality MapReject initial interference phase gradient estimation figure In the lower interferometric phase gradient estimated value of quality, obtain final interferometric phase gradient estimation figure
WithMethod is:
3.1, it will estimate Quality MapMiddle all elements are normalized, i.e.,
Wherein max () expression is maximized operation；
3.2, quality threshold is setFor normalizing Quality MapIn be higher than quality threshold element,
Then determine initial interference phase gradient estimation figureWithIt is credible, final interferometric phase gradient estimation figure
3.3, for normalizing Quality MapIn each be lower than quality threshold element, it is believed that initial interference phase
Gradient estimation figureWithIt is insincere, it needs to reevaluate interferometric phase gradient, finally interferes phase
Potential gradient estimation figure are as follows:
4th step utilizes obtained final interferometric phase gradient estimation figureWithCompensation filter window
Interferometric phase gradient in mouthful, to initial interferometric phase image s again1(x, y) carries out mean filter, multiple interferometric phase after being filtered
Scheme s3(x, y), method are: initially answering interferometric phase image s1Filter window, such as window size 7 × 7, filtering are set in (x, y)
Any pixel s in interferometric phase image is answered afterwards3(x, y) are as follows:
Wherein, N is sum of all pixels in filter window, and x and y are respectively the two dimension of any pixel in initial interferometric phase image again
Coordinate, fx' it is initially to answer pixel s in interferometric phase image1The final interferometric phase gradient estimated value of (x, y) corresponding X direction,
fy' it is initially to answer pixel s in interferometric phase image1The final interferometric phase gradient estimated value of (x, y) corresponding y direction.
Embodiment
The present embodiment uses the Italian volcano the Etna region InSAR data enrolled the SIR-C/X-SAR month, main SAR data
Magnitude image it is as shown in Figure 4.
The first step carries out interference processing to the major-minor SAR complex pattern after accuracy registration, generates initially interferometric phase image s again1
(x, y) is as shown in figure 5, to interferometric phase image s is initially answered1(x, y) carries out 3 × 3 mean filters, and obtained coarse filtration wave interferes phase again
Bitmap s2(x, y) is as shown in fig. 6, by 3 × 3 mean filters, while retaining original interference phase diagram detailed information, substantially
The phase noise in original interference phase diagram is inhibited, the robustness of subsequent intervention phase gradient estimation is improved.
Second step answers interferometric phase image s in coarse filtration wave2Setting interferometric phase gradient estimates window in (x, y), is based on FFT
Estimate two-dimentional initial interference phase gradient estimation figureRespectively as shown in Figure 7 and Figure 8, estimated
It is as shown in Figure 9 to count Quality Map.
Third step, according to estimation Quality MapReject initial interference phase gradient estimation figure In the lower interferometric phase gradient estimated value of quality, obtain final interferometric phase gradient estimation figureWithRespectively as shown in Figure 10 and Figure 11, estimate that the quality in figure is lower dry by rejecting initial interference phase gradient
Phase gradient estimated value is related to, interferometric phase gradient estimated accuracy is improved.
4th step utilizes obtained final interferometric phase gradient estimation figureWithCompensation filter window
In interferometric phase gradient, to initially again interferometric phase image s1(x, y) carries out 7 × 7 mean filters, interferes phase after being filtered again
Bitmap s3(x, y), as shown in figure 12.Comparison diagram 1 and Figure 12 are it can be seen that the method for the present invention is keeping interferometric phase detailed information
While, effectively reduce Random phase noise.By statistics, initially residual points are 202002 in interferometric phase image again,
Residual points are 4712 in multiple interferometric phase image after filtering, and residual rejection rate reaches 97.6%, and residual point refers to for interferometric phase here
Four adjacent pixels, the closed loop integral for calculating its interferometric phase gradient are not zero in figure.
The present embodiment is first to initially interferometric phase image carries out 3 × 3 mean filters again, to improve the estimation of interferometric phase gradient
Robustness；When obtaining initial interference phase gradient estimation figure, while obtaining the estimation Quality Map of interferometric phase gradient, quality
Coefficient is higher than the interferometric phase gradient initial value of threshold value as final estimated value, and quality is picked lower than the initial value of threshold value
It removes, the final estimated value of interferometric phase gradient of the position is weighted using the higher interferometric phase gradient initial value of ambient quality
Averagely obtain.Finally, using the interferometric phase gradient in the interferometric phase gradient compensating interferometer phase filtering window for estimating to obtain,
Mean filter is carried out again, completes the whole process of InSAR phase filtering.
The present embodiment additionally provides a kind of InSAR filtering system that interferometric phase gradient is estimated based on FFT high-precision, packet
Include: the first module generates initially interferometric phase image s again for carrying out interference processing to the major-minor SAR complex pattern after registration1(x,
Y), then to s1(x, y) progress mean filter obtains coarse filtration wave and answers interferometric phase image s2(x,y)；Wherein, x and y is interferometric phase
The two-dimensional coordinate of any pixel in figure；Second module, for answering interferometric phase image s in coarse filtration wave2Interferometric phase is set in (x, y)
Gradient estimates window, estimates to obtain initial interference phase gradient estimation figure (two-dimentional interferometric phase gradient initial value) based on FFTWithEstimate Quality Map (quality coefficient)Wherein, p and q is that initial interference phase gradient is estimated
The two-dimensional coordinate of any pixel in figure is counted,Estimate to scheme for X direction initial interference phase gradient,For
Y direction initial interference phase gradient estimation figure；Third module, for according to estimation Quality Map (quality coefficient)
Reject initial interference phase gradient estimation figureIn quality lower interferometric phase gradient estimation
Value, obtains final interferometric phase gradient estimation figureWith4th module, for utilizing final interferometric phase
Gradient estimation figureWithInterferometric phase gradient in compensation filter window, to initial interferometric phase image again
s1(x, y) carries out mean filter, multiple interferometric phase image s after being filtered3(x,y)。
The present embodiment solves the contradiction between the filtering view number faced in InSAR phase filtering and terrain information holding,
Overcome coherence it is poor when the interferometric phase gradient estimation difficult points such as robustness difference, provide branch for the high-precision InSAR measurement of higher degree
Support has important application value.
Embodiment described above is the present invention more preferably specific embodiment, and those skilled in the art is in this hair
The usual variations and alternatives carried out in bright technical proposal scope should be all included within the scope of the present invention.
Claims (10)
1. a kind of InSAR filtering method for estimating interferometric phase gradient based on FFT high-precision, which is characterized in that the method packet
Include following steps:
Step (1): carrying out interference processing to the major-minor SAR complex pattern after registration, generates initially interferometric phase image s again1(x, y), so
Afterwards to s13 × 3 mean filters of (x, y) progress obtain coarse filtration wave and answer interferometric phase image s2(x,y)；Wherein, x and y is respectively initial multiple
The two-dimensional coordinate of any pixel in interferometric phase image；
Step (2): interferometric phase image s is answered in coarse filtration wave2Setting interferometric phase gradient estimates window in (x, y), is estimated based on FFT
Obtain initial interference phase gradient estimation figureWithEstimate Quality MapWherein, p and q is first
Beginning interferometric phase gradient estimates the two-dimensional coordinate of any pixel in figure,Estimate for X direction initial interference phase gradient
Meter figure,Estimate to scheme for y direction initial interference phase gradient；
Step (3): according to estimation Quality MapReject initial interference phase gradient estimation figure
In quality be lower than 0.4 interferometric phase gradient estimated value, obtain final interferometric phase gradient estimation figureWith
Step (4): final interferometric phase gradient estimation figure is utilizedWithInterference in compensation filter window
Phase gradient, to initial interferometric phase image s again1(x, y) carries out mean filter, multiple interferometric phase image s after being filtered3(x,y)。
2. the InSAR filtering method according to claim 1 for estimating interferometric phase gradient based on FFT high-precision, feature
It is: in the step (1), to initial interferometric phase image s again1The method that (x, y) carries out 3 × 3 mean filters are as follows: initial
Filter window is set in multiple interferometric phase image, and window size is 3 × 3, and it is average then to carry out plural number to pixel in window.
3. the InSAR filtering method according to claim 2 for estimating interferometric phase gradient based on FFT high-precision, feature
Be: coarse filtration wave answers interferometric phase image s2(x, y) isN is spectral window
Sum of all pixels in mouthful, -1≤lx,ly≤ 1 is two integer variables, with lxAnd lyValue is different, (x+lx,y+ly) indicate with (x,
It y) is the different location in the filter window of center, x and y are respectively the initial two-dimensional coordinate for answering any pixel in interferometric phase image.
4. the InSAR filtering method according to claim 1 for estimating interferometric phase gradient based on FFT high-precision, feature
Be: the step (2) specifically comprises the following steps:
(2.1) interferometric phase image s is answered in coarse filtration wave2Setting interferometric phase gradient estimates window, window size N in (x, y)f×
Nf, window center is located at (x, y), and the coarse filtration wave in window interferes diagram data to be g (m, n) again；
(2.2) zero padding is first carried out to g (m, n) to operate to obtain gz(m, n), then to gz(m, n) carries out two-dimensional fast fourier transform
Obtain 2-d spectrum Gz(u,v)；
(2.3) 2-d spectrum GzProminent frequency coordinate is (u in (u, v)max,vmax), calculate fx、fyAnd C:
fx=2 π (umax-Nz/2)/Nz
fy=2 π (vmax-Nz/2)/Nz
Wherein, | | it is operated for calculated complex modulus value, NzThe g operated for zero paddingzThe two-dimensional of (m, n), i.e. gz(m,n)
Data scale be Nz×Nz, Gz(umax,vmax) it is to gz(m, n) carries out two-dimensional fast fourier transform and obtains 2-d spectrum Gz
The frequency spectrum of power maximum, u in (u, v)maxAnd vmaxRespectively 2-d spectrum GzThe two-dimensional frequency of power maximum in (u, v)
Value, fxAnd fyThe X direction and y direction initial interference phase of all pixels respectively in interferometric phase gradient estimation window
Gradient estimated value, C are that interferometric phase gradient estimates that the interferometric phase gradient of all pixels in window estimates quality coefficient；
(2.4) position for changing interferometric phase gradient estimation window, so that new window center is located at (x+Nf, y) or (x, y+
Nf), it repeats step (2.1) and arrives (2.3), until completing the interferometric phase ladder that estimation coarse filtration wave answers each pixel in interferometric phase image
Degree and quality coefficient, obtain initial interference phase gradient estimation figureWithAnd estimation Quality Map
5. the InSAR filtering method according to claim 1 for estimating interferometric phase gradient based on FFT high-precision, feature
Be: the step (3) specifically comprises the following steps:
(3.1) it will estimate Quality MapMiddle all elements are normalized；
(3.2) quality threshold is setFor normalizing Quality MapIn be higher than quality threshold element, then
Determine initial interference phase gradient estimation figureWithEffectively, final interferometric phase gradient estimation figure
(3.3) for normalizing Quality MapIn each be lower than quality threshold element, then determine initial interference phase
Gradient estimation figureWithIn vain, interferometric phase gradient is reevaluated, final interferometric phase gradient
Estimation figure are as follows:
Wherein,The weighting weight of corresponding element in figure, k are estimated for initial interference phase gradientpAnd kqFor integer, indicate most
Two-dimensional coordinate when whole interferometric phase gradient is estimated relative to (p, q),For (p+ in estimation Quality Map
kp,q+kq) at estimation quality coefficient.
6. the InSAR filtering method according to claim 1 for estimating interferometric phase gradient based on FFT high-precision, feature
It is: in the step (4), multiple interferometric phase image s after filtering3(x, y) are as follows:
Wherein, it is sum of all pixels in filter window that window size, which is 7 × 7, N, and x and y are respectively initially to answer in interferometric phase image to appoint
The two-dimensional coordinate for pixel of anticipating, f 'xInitially to answer pixel s in interferometric phase image1The final interferometric phase of (x, y) corresponding X direction
Gradient estimated value, f 'yInitially to answer pixel s in interferometric phase image1The final interferometric phase gradient of (x, y) corresponding y direction is estimated
Evaluation.
7. the InSAR filtering method according to claim 4 for estimating interferometric phase gradient based on FFT high-precision, feature
It is: in step (2.1), interferes the formula of diagram data g (m, n) as follows again in window:
G (m, n)=s2(x+m-Nf/2,y+n-Nf/ 2),
Wherein, 1≤m, n≤NfFor positive integer, m and n indicate that interferometric phase gradient estimates the two-dimensional coordinate of pixel in window.
8. the InSAR filtering method according to claim 4 for estimating interferometric phase gradient based on FFT high-precision, feature
It is: in step (2.2), the g after zero paddingz(m, n) points are Nz×Nz, Nz=128, gz(m, n) is indicated are as follows:
9. the InSAR filtering method according to claim 5 for estimating interferometric phase gradient based on FFT high-precision, feature
It is: in step (3.1), will estimates Quality MapMiddle all elements, which are normalized, to be obtained by the following formula:
Wherein, max () expression is maximized operation.
10. a kind of InSAR filtering system for estimating interferometric phase gradient based on FFT high-precision, characterized by comprising:
First module generates initially interferometric phase image s again for carrying out interference processing to the major-minor SAR complex pattern after registration1(x,
Y), then to s1(x, y) progress mean filter obtains coarse filtration wave and answers interferometric phase image s2(x,y)；Wherein, x and y is initial multiple dry
Relate to the two-dimensional coordinate of any pixel in phase diagram；
Second module, for answering interferometric phase image s in coarse filtration wave2Setting interferometric phase gradient estimates window in (x, y), is based on FFT
Estimation obtains initial interference phase gradient estimation figure (two-dimentional interferometric phase gradient initial value)WithEstimate
It counts Quality Map (quality coefficient)Wherein, p and q is that initial interference phase gradient estimates that the two dimension of any pixel in figure is sat
Mark,Estimate to scheme for X direction initial interference phase gradient,For y direction initial interference phase ladder
Degree estimation figure；
Third module, for according to estimation Quality Map (quality coefficient)Reject initial interference phase gradient estimation figureIn the lower interferometric phase gradient estimated value of quality, obtain the estimation of final interferometric phase gradient
FigureWith
4th module, for utilizing final interferometric phase gradient estimation figureWithIn compensation filter window
Interferometric phase gradient, to initial interferometric phase image s again1(x, y) carries out mean filter, multiple interferometric phase image s after being filtered3(x,y)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810779761.XA CN109143234B (en) | 2018-07-16 | 2018-07-16 | InSAR filtering method and system for estimating interference phase gradient based on FFT high precision |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810779761.XA CN109143234B (en) | 2018-07-16 | 2018-07-16 | InSAR filtering method and system for estimating interference phase gradient based on FFT high precision |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109143234A true CN109143234A (en) | 2019-01-04 |
| CN109143234B CN109143234B (en) | 2021-03-26 |
Family
ID=64800679
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810779761.XA Active CN109143234B (en) | 2018-07-16 | 2018-07-16 | InSAR filtering method and system for estimating interference phase gradient based on FFT high precision |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109143234B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117688341A (en) * | 2024-01-31 | 2024-03-12 | 安徽水安建设集团股份有限公司 | Deep foundation pit detection system and method based on BIM technology |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2942637A1 (en) * | 2014-05-09 | 2015-11-11 | NEC Corporation | Change detection device, change detection method and recording medium |
| CN105116410A (en) * | 2015-07-20 | 2015-12-02 | 西北农林科技大学 | Interferometric phase adaptive filtering algorithm based on linear model matching |
| CN105929399A (en) * | 2016-04-25 | 2016-09-07 | 电子科技大学 | Interference SAR data imaging and elevation estimation method |
| US10006995B2 (en) * | 2014-08-04 | 2018-06-26 | University Of Seoul Industry Cooperation Foundation | Method and apparatus for stacking multi-temporal MAI interferograms |
-
2018
- 2018-07-16 CN CN201810779761.XA patent/CN109143234B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2942637A1 (en) * | 2014-05-09 | 2015-11-11 | NEC Corporation | Change detection device, change detection method and recording medium |
| US10006995B2 (en) * | 2014-08-04 | 2018-06-26 | University Of Seoul Industry Cooperation Foundation | Method and apparatus for stacking multi-temporal MAI interferograms |
| CN105116410A (en) * | 2015-07-20 | 2015-12-02 | 西北农林科技大学 | Interferometric phase adaptive filtering algorithm based on linear model matching |
| CN105929399A (en) * | 2016-04-25 | 2016-09-07 | 电子科技大学 | Interference SAR data imaging and elevation estimation method |
Non-Patent Citations (2)
| Title |
|---|
| BINGQIAN CHEN: "Combining SAR interferometric phase and intensityinformation for monitoring of large gradient deformation in coal mining area", 《EUROPEAN JOURNAL OF REMOTE SENSING》 * |
| 蒋锐: "一种基于等效残差点的InSAR相位解缠绕方法", 《南京航空航天大学学报》 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117688341A (en) * | 2024-01-31 | 2024-03-12 | 安徽水安建设集团股份有限公司 | Deep foundation pit detection system and method based on BIM technology |
| CN117688341B (en) * | 2024-01-31 | 2024-05-14 | 安徽水安建设集团股份有限公司 | Deep foundation pit detection system and method based on BIM technology |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109143234B (en) | 2021-03-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105528785B (en) | A kind of binocular vision image solid matching method | |
| CN106373087B (en) | A kind of image super-resolution rebuilding method improving initial estimation | |
| Dai et al. | An accurate phase unwrapping algorithm based on reliability sorting and residue mask | |
| KR101700928B1 (en) | Bayer pattern image demosaicking method and apparatus based on multi-directional weighted interpolation and guided filter | |
| CN109633648A (en) | A kind of more baseline phase estimation devices and method based on possibility predication | |
| AU2005280315A1 (en) | Method and system for motion correction in a sequence of images | |
| CN101419669B (en) | Three-dimensional human ear extracting method based on profile wave convert | |
| CN109064418A (en) | A kind of Images Corrupted by Non-uniform Noise denoising method based on non-local mean | |
| CN107014313A (en) | The method and system of weighted least-squares phase unwrapping based on S-transformation ridge value | |
| CN109143234A (en) | The InSAR filtering method and system of interferometric phase gradient are estimated based on FFT high-precision | |
| CN107290720A (en) | A kind of detection method of the moving spot targets based on passive microwave interference imaging | |
| CN109443250B (en) | Structured light three-dimensional surface shape vertical measurement method based on S transformation | |
| CA2423377C (en) | Noise reduction in images | |
| CN113269691B (en) | SAR image denoising method for noise affine fitting based on convolution sparsity | |
| CN111179333A (en) | Defocus fuzzy kernel estimation method based on binocular stereo vision | |
| WO2013134998A1 (en) | 3d reconstruction method and system | |
| Gilman et al. | Least-squares optimal interpolation for fast image super-resolution | |
| Kelkar et al. | AmbientFlow: Invertible generative models from incomplete, noisy measurements | |
| CN105116410B (en) | The interferometric phase image adaptive filter algorithm matched based on linear model | |
| CN106707283B (en) | Phase-unwrapping algorithm based on tasteless information filter | |
| Islam et al. | Super resolution of 3d MRI images using a Gaussian scale mixture model constraint | |
| CN110298870A (en) | Processing method, processing unit and the terminal of image | |
| Gol Gungor et al. | A subspace‐based coil combination method for phased‐array magnetic resonance imaging | |
| CN110440935B (en) | Phase unwrapping method based on extended information filtering | |
| Ferrari et al. | Distributed image reconstruction for very large arrays in radio astronomy |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |