Ground-based interferometric synthetic aperture radar-based atmospheric phase compensation method
The present invention relates to the technical field of synthetic aperture radar, and in particular to a ground-based interference synthetic aperture radar-based atmospheric phase compensation method. The present invention first uses the inverse fast Fourier transform algorithm to rapidly and effectively realize the focusing in the range and cross-range dimension. Then a triple threshold method which combines coherence coefficient, amplitude, and amplitude dispersion index is used to select reliable PS points. Finally, under a full consideration of the spatial correlation of the atmospheric phase, the atmospheric phase is estimated by using a two-dimensional-polynomial model. The present invention can rapidly and accurately estimate and compensate the atmospheric phase, is helpful in improving the accuracy of GB-InSAR real-time measurement, and valuable and universal in practical application.
1. A ground-based interferometric synthetic aperture radar (GB-InSAR)-based atmospheric phase compensation method, comprising:
a. imaging processing: setting as using a step frequency continuous wave (SFCW) technique, transforming a range compression from a frequency domain to a spatial domain by an inverse fast Fourier transform, regarding a sampling of different observation points in ground-based synthetic aperture radar (GB-SAR) along a linear rail as a frequency domain sampling with a linear frequency shift, performing cross-range focusing by the inverse fast Fourier transform as well, and adding windows to an image after focusing to obtain a final two-dimensional synthetic aperture radar (SAR) image;
b. image registration: firstly, matching the final two-dimensional SAR images obtained by continuous observations; sequentially, finding out a transformation from a slave image S geometry to a master image M geometry; finally, resampling a slave image S to a master image M geometry, wherein step b specifically comprises:
b1. image matching: selecting a series of pixel points (x,y) from the master image M, and estimating a shift k xy =(k x ,k y ) of the image in a range direction and a cross-range direction according to a cross-correlation function C(k xy ) of a square of amplitude values of the master image and the slave image, wherein the cross-correlation function C(k xy ) is expressed as follows:
C
(
k
xy
)
=
∑
i
,
j
∈
H
xy
[
M
(
i
,
j
)
-
μ
M
]
[
S
(
i
-
k
x
,
j
-
k
y
)
-
μ
S
]
σ
M
σ
S
(
formula
1
)
in formula 1, H xy represents a window function, μ M , μ S , σ M , σ S represent mean values and variances of images in the window, respectively; the maximum value of the cross-correlation function C(k xy ) is used to represent the shift k xy =(k x ,k y ) between the master image and the slave image;
b2. transformation estimation: estimating transformation parameters p x (x,y),p y (x,y) between the master image and the slave image by the following two-dimensional-polynomial:
{
p
x
(
x
,
y
)
=
a
x
+
b
x
x
+
c
x
y
+
d
x
x
2
+
e
x
xy
+
f
x
y
2
p
y
(
x
,
y
)
=
a
y
+
b
y
x
+
c
y
y
+
d
y
x
2
+
e
y
xy
+
f
y
y
2
(
formula
2
)
wherein a x , b x , c x , d x , e x , f x and a y , b y , c y , d y , e y , f y are estimated by the least squares method, and an observation model is expressed as follows:
{
x
S
=
x
M
+
p
x
(
x
,
y
)
+
ξ
x
y
S
=
y
M
+
p
y
(
x
,
y
)
+
ξ
y
(
formula
3
)
in formula 3, (x S ,y S ) is a coordinate on the slave image corresponding to the point (x M ,y M ) on the master image, and (ξ x ,ξ y ) is a residual observation error; a representation of the shift (k x ,k y ) obtained is expressed as follows:
{
k
x
=
p
x
(
x
,
y
)
+
ξ
x
k
y
=
p
y
(
x
,
y
)
+
ξ
y
(
formula
4
)
b3. image resampling: using a truncated Sinc as an interpolation signal, transforming the slave image to the master image, namely, each pixel point of the master image and the slave image matches a same portion of a terrain of a measured area;
c. interferometric phase filtering and phase unwrapping: performing a differential interference processing on the master image and the slave image registered in step b to obtain an interference image; for the GB-SAR system, an interference phase of an α th pixel point in an interference image is expressed as follows:
Δψ a =Δψ def,a +Δψ atm,a +Δψ noise,a +2 nπ (formula 5)
wherein in formula 5, Δψ def,a represents a phase caused by a deformation of target in two observations, Δψ atm,a represents a phase difference caused by atmospheric factors in the two observations, and Δψ noise,a represents a phase difference caused by thermal noise factors;
a minimum cost flow algorithm is used for phase unwrapping, and a periodic mean filtering method is used to filter out an interference of noise before the phase unwrapping;
d. permanent scatterer (PS) point selection: grouping SAR images, wherein every F SAR images are set as a group and F denotes the number of SAR images; firstly, the PS candidate selection is performed by using amplitude and amplitude dispersion information of each group; secondly, a secondary screening is performed by using correlation coefficient information of all PS candidates and selecting out a final PS point set step d further comprises:
d1. PS candidate selection: selecting PS candidates by using an amplitude and an amplitude dispersion index;
d11. selecting the PS candidates by using amplitude dispersion information, wherein a calculation formula of an amplitude dispersion is as follows:
Υ
amp
=
σ
amp
m
amp
(
formula
6
)
In formula 6, σ amp represents amplitude standard deviation of pixel point in F images, and m amp represents amplitude average of the pixel point in F images; in a process of an amplitude dispersion selection, an amplitude dispersion threshold ρ adi is set firstly, then pixel points having the amplitude dispersion less than the threshold ρ adi are selected as the PS candidates;
d12. further selecting the PS candidates by using the amplitude information, wherein an amplitude selection method is expressed by a mathematical formula as follows:
{
β
≥
β
thr
,
PS
candidate
β
<
β
thr
,
reject
point
(
formula
7
)
in formula 7, β is average amplitude information of the PS candidates selected according to the amplitude dispersion information in F images participating in the selection, and β thr is an amplitude threshold;
d2. secondary screening: selecting a final PS point set by using coherence coefficient information of all PS candidates, and a mathematical expression of the coherence coefficient is as follows:
γ
=
∑
ξ
∑
ζ
M
·
S
*
∑
ξ
∑
ζ
M
2
∑
ξ
∑
ζ
S
2
(
formula
8
)
in formula 8, |·| and (·)* respectively represent modulus operation and conjugate operation, and ξ×ζ denotes a size of a two-dimensional window function; pixel points having the coherence coefficient less than a threshold ρ coh are removed by setting the coherence coefficient threshold ρ coh , and the remaining PS candidates are final selection results of PS points;
e. atmospheric phase compensation: estimating the atmospheric phase by using a two-dimensional-polynomial model based on the PS points selected in step d, and establishing a two-dimensional-polynomial model with a degree transformation of r:
Δψ
atm
(
p
,
q
)
=
b
00
+
∑
t
=
0
r
b
t
,
(
r
-
t
)
·
(
p
-
p
ref
)
t
·
(
q
-
q
ref
)
r
-
t
(
formula
9
)
wherein in formula 9, represents an atmospheric phase, (p,q) represents a coordinate of a pixel point, (p ref ,q ref ) represents a coordinate of a reference pixel point, and b t,(r-t) represents a polynomial model parameter; an atmospheric phase screen (APS) model parameter is estimated based on the PS points selected, a deformation phase of the PS point selected is set as zero, so the following equation can be obtained:
Δψ( p,q )=Δψ atm ( p,q )+Δψ noise ( p,q )
an estimated model parameter {circumflex over (b)} t,(r-t) is estimated by a least squares method, then the atmospheric phase of a deformation area is calculated by using the estimated model parameter {circumflex over (b)} t,(r-t) ; after the atmospheric phase compensation, a true phase value caused by deformation is obtained, so that a deformation of a monitoring area is calculated.
2. The GB-InSAR-based atmospheric phase compensation method according to claim 1 , wherein in step d11, the PS candidates are selected by using the amplitude dispersion information, the amplitude dispersion threshold ρ adi is set as 0.1.
3. The GB-InSAR-based atmospheric phase compensation method according to claim 1 , wherein in step d2, the coherence coefficient threshold ρ coh is set as 0.9 during the secondary screening.