US 10,705,205 B2 · Ground-based interferometric synthetic aperture radar-based atmospheric phase compensation method

Abstract

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.

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 kxy=(kx,ky) of the image in a range direction and a cross-range direction according to a cross-correlation function C(kxy) of a square of amplitude values of the master image and the slave image, wherein the cross-correlation function C(kxy) 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, Hxy 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(kxy) is used to represent the shift kxy=(kx,ky) between the master image and the slave image;
b2. transformation estimation:
estimating transformation parameters px (x,y),py (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 ax, bx, cx, dx, ex, fx and ay, by, cy, dy, ey, fy 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, (xS,yS) is a coordinate on the slave image corresponding to the point (xM,yM) on the master image, and (ξx,ξy) is a residual observation error;
a representation of the shift (kx,ky) 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+2nπ (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 mamp 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, (pref,qref) represents a coordinate of a reference pixel point, and bt,(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.

385 more elements

Ground-based interferometric synthetic aperture radar-based atmospheric phase compensation method

Life of this patent

Ground-based interferometric synthetic aperture radar-based atmospheric phase compensation method

Life of this patent

Abstract

Claims3