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.

Claims

1. A GB-InSAR-based atmospheric phase compensation method, comprising: a. imaging processing: setting as using an 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 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 SAR image; b. image registration: registering the final two-dimensional SAR images obtained by continuous observations, sequentially, namely, finding out a transformation from the slave image S geometry to the master image M geometry, and resampling the slave image S to the 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.sub.xy=(k.sub.x,k.sub.y) of the image in a range direction and a cross-range direction according to a cross-correlation function C(k.sub.xy) of a square of amplitude values of the master image and the slave image, wherein the cross-correlation function C(k.sub.xy) is expressed as follows: $\begin{array}{cc}C\left(k_{\mathrm{xy}}\right)=\frac{\underset{i,j{H}_{\mathrm{xy}}}{.Math.}.Math.\left[M\left(i,j\right)-{}_M\right]\left[S\left(i-k_x,j-k_y\right)-{}_S\right]}{{}_M.Math.{}_S}& \left(\mathrm{formula}.Math..Math.1\right)\end{array}$ in formula 1, H.sub.xy represents a window function, .sub.M, .sub.S, .sub.M, .sub.S represent mean values and variances of images in the window, respectively; the maximum value of the cross-correlation function C(k.sub.xy) is used to represent the shift k.sub.xy=(k.sub.x,k.sub.y) between the master image and the slave image; b2. transformation estimation: estimating transformation parameters p.sub.x(x,y),p.sub.y(x,y) between the master image and the slave image by the following two-dimensional-polynomial: $\begin{array}{cc}\{\begin{array}{c}{p}_x\left(x,y\right)=a_x+b_x.Math.x+c_x.Math.y+d_x.Math.x^2+e_x.Math.\mathrm{xy}+f_x.Math.y^2\\ p_y\left(x,y\right)=a_y+b_y.Math.x+c_y.Math.y+d_y.Math.x^2+e_y.Math.\mathrm{xy}+f_y.Math.y^2\end{array}& \left(\mathrm{formula}.Math..Math.2\right)\end{array}$ wherein a.sub.x, b.sub.x, c.sub.x, d.sub.x, e.sub.x, f.sub.x and a.sub.y, b.sub.y, c.sub.y, d.sub.y, e.sub.y, f.sub.y are estimated by the least squares method, and an observation model is expressed as follows: $\begin{array}{cc}\{\begin{array}{c}{x}_S=x_M+p_x\left(x,y\right)+{}_x\\ y_S=y_M+p_y\left(x,y\right)+{}_y\end{array}& \left(\mathrm{formula}.Math..Math.3\right)\end{array}$ in formula 3, (x.sub.S,y.sub.S) is a coordinate on the slave image corresponding to the point (x.sub.M,y.sub.M) on the master image, and (.sub.x,.sub.y) is a residual observation error; a representation of the shift (k.sub.x,k.sub.y) obtained is expressed as follows: $\begin{array}{cc}\{\begin{array}{c}{k}_x=p_x\left(x,y\right)+{}_x\\ k_y=p_y\left(x,y\right)+{}_y\end{array}& \left(\mathrm{formula}.Math..Math.4\right)\end{array}$ 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 GB-SAR system, an interference phase of an .sup.th pixel point in an interference image is expressed as follows:
.sub.a=.sub.def,a+.sub.atm,a+.sub.noise,a+2n(formula 5) wherein in formula 5, .sub.def,a represents a phase caused by a deformation of target in two observations, .sub.atm,a represents a phase difference caused by atmospheric factors in the two observations, and .sub.noise,a represents a phase difference caused by factors such as thermal noise; 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. PS point selection: grouping SAR images, wherein every F SAR images are set as a group, and the PS candidate selection is performed for each group; using correlation coefficient information of all PS candidates to perform a secondary screening for 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: $\begin{array}{cc}{}_{\mathrm{amp}}=\frac{{}_{\mathrm{amp}}}{m_{\mathrm{amp}}}& \left(\mathrm{formula}.Math..Math.6\right)\end{array}$ In formula 6, .sub.amp represents a time sequential amplitude standard deviation of target point, and m.sub.amp represents a time sequential amplitude average of the target point; in a process of an amplitude dispersion selection, an amplitude dispersion threshold .sub.adi is set first, then pixel points having the amplitude dispersion less than the threshold .sub.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: $\begin{array}{cc}\{\begin{array}{c}{}_{\mathrm{thr}},\mathrm{PS}.Math..Math.\mathrm{candidate}\\ <{}_{\mathrm{thr}},\mathrm{reject}.Math..Math.\mathrm{point}\end{array}& \left(\mathrm{formula}.Math..Math.7\right)\end{array}$ 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 .sub.thr is an amplitude threshold; d2. performing a secondary screening by using coherence coefficient information, wherein a window with a size of is used to calculate a coherence coefficient, and a mathematical expression of the coherence coefficient is as follows: $\begin{array}{cc}=\frac{\underset{}{.Math.}.Math.\underset{}{.Math.}.Math.M.Math.S^{*}}{\sqrt{\underset{}{.Math.}.Math.\underset{}{.Math.}.Math.{.Math.M.Math.}^2.Math.\underset{}{.Math.}.Math.\underset{}{.Math.}.Math.{.Math.S.Math.}^2}}& \left(\mathrm{formula}.Math..Math.8\right)\end{array}$ the coherence coefficient of each of the PS candidates in a current image and a previous image is calculated, then all points having the coherence coefficient less than a threshold .sub.coh are removed by setting the coherence coefficient threshold .sub.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: $\begin{array}{cc}{}_{\mathrm{atm}}\left(p,q\right)=b_{00}+\underset{t=0}{\overset{r}{.Math.}}.Math.b_{t,\left(r-t\right)}.Math.{\left(p-p_{\mathrm{ref}}\right)}^t.Math.{\left(q-q_{\mathrm{ref}}\right)}^{r-t}& \left(\mathrm{formula}.Math..Math.9\right)\end{array}$ wherein in formula 9, .sub.atmo represents an atmospheric phase, (p, q) represents a coordinate of a pixel point, (p.sub.ref,q.sub.ref) represents a coordinate of a reference pixel point, and b.sub.t,(r-t) represents a polynomial model parameter; an 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)=.sub.atm(p,q)+.sub.noise(p,q) an estimated model parameter {circumflex over (b)}.sub.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)}.sub.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 .sub.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 .sub.coh is set as 0.9 during the secondary screening.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0030] FIG. 1 is a flow chart of the present invention.

[0031] FIG. 2 shows the result of PS point selection of the data measured in Yanguo gorge, Yongjing County, Gansu province of China.

[0032] FIG. 3 shows the result of the data measured in Yanguo gorge Yongjing County, Gansu Province of China before and after the atmospheric phase compensation.

[0033] FIG. 4 shows the unbiased estimation result of the atmospheric phase variance error of the data measured in Yanguo gorge, Yongjing County, Gansu Province of China.

DETAILED DESCRIPTION

[0034] The present invention will be further described below with reference to the drawings of the specification.

[0035] The processing flow chart of the GB-InSAR-based atmospheric phase compensation method is shown in FIG. 1. The processing flow consists of five steps: imaging processing, image registration, interferometric phase filtering and phase unwrapping, PS point selection, and atmospheric phase compensation, which are specifically described as follows.

[0036] Step(1), imaging processing: because SFCW technique is used, the range compression can be transformed from the frequency domain to the spatial domain by the inverse fast Fourier transform. From the other aspect, since the sampling of different observation points in GB-SAR along the linear rail can be regarded as the frequency domain sampling with a linear frequency shift, the cross-range focusing can be realized by the inverse fast Fourier transform as well. To reduce the side lobes, a window function can be applied after focusing to obtain the final two-dimensional SAR image. The optimal window function used in the present invention is the Hamming window:

W.sub.Hamming(l,a)=(1) cos[2l(L1)],0lL1

[0037] L is the length of the window, and generally, the value of a is 0.46.

[0038] Step (2), image registration: the two-dimensional SAR images obtained by continuous observations are sequentially registered. To perform the image registration, the transformation from the slave image S geometry to the master image M geometry has to be found and then the slave image S has to be resampled to the master image M geometry, which specifically includes the following steps.

[0039] 2.1) Image matching: a series of pixel points (x,y) are selected from the master image M, and the shift k.sub.xy=(k.sub.x,k.sub.y) of the image in range and cross-range direction is estimated according to the cross-correlation function C(k.sub.xy) of the square of the amplitude values of the two images. The cross-correlation function C(k.sub.xy) is expressed as follows:

[00009]$C\left(k_{\mathrm{xy}}\right)=\frac{\underset{i,j{H}_{\mathrm{xy}}}{.Math.}.Math.\left[M\left(i,j\right)-{}_M\right]\left[S\left(i-k_x,j-k_y\right)-{}_S\right]}{{}_M.Math.{}_S}$

[0040] where H.sub.xy represents a window function, and .sub.M, .sub.S, .sub.M, .sub.S as represent the mean values and variances of images in the window, respectively. The maximum value of cross-correlation function C(k.sub.xy) is used to represent the shift k.sub.xy=(k.sub.x,k.sub.y) between two images.

[0041] 2.2) Transformation estimation: the transformation parameters p.sub.x(x, y), p.sub.y(x, y) between two images are estimated by the following two dimensional-polynomial:

[00010]$\{\begin{array}{c}{p}_x\left(x,y\right)=a_x+b_x.Math.x+c_x.Math.y+d_x.Math.x^2+e_x.Math.\mathrm{xy}+f_x.Math.y^2\\ p_y\left(x,y\right)=a_y+b_y.Math.x+c_y.Math.y+d_y.Math.x^2+e_y.Math.\mathrm{xy}+f_y.Math.y^2\end{array}$

[0042] where, a.sub.x, b.sub.x, c.sub.x, d.sub.x, e.sub.x, f.sub.x and a.sub.y, b.sub.y, c.sub.y, d.sub.y, e.sub.y, f.sub.y are estimated by the least squares method. The observation model is:

[00011]$\{\begin{array}{c}{x}_S=x_M+p_x\left(x,y\right)+{}_x\\ y_S=y_M+p_y\left(x,y\right)+{}_y\end{array}$

[0043] where, (x.sub.s,y.sub.s) is a coordinate on the slave image corresponding to the point (x.sub.M, y.sub.M) on the master image, and (.sub.x, .sub.y) is the residual observation error. The representation of the shift (k.sub.x, k.sub.y) obtained is expressed as follows:

[00012]$\{\begin{array}{c}{k}_x=p_x\left(x,y\right)+{}_x\\ k_y=p_y\left(x,y\right)+{}_y\end{array}$

[0044] 2.3) Image resampling: in the present invention, the truncated Sinc interpolation is used, and the truncated Sinc kernel size is the key parameter. The typical interpolation kernel size is 8 pixels. The output of this step is the transformation from the slave image to the master image, namely, each pixel point of the master and slave images matches the same portion of terrain of the measured area.

[0045] Step (3), interferometric phase filtering and phase unwrapping: the differential interference processing is performed on the two images registered in step (2) to obtain an interference image. For GB-SAR system, the interference phase of the .sup.th pixel point in the interference image can be expressed as follows:

.sub.a=.sub.def,a+.sub.atm,a+.sub.noise,a+2n

[0046] where, .sub.def,a represents the phase caused by the deformation of target in two observations, .sub.atm,a represents the phase difference caused by atmospheric factors in two observations, and .sub.noise,a represents the phase difference caused by factors such as thermal noise etc. As a result of the periodicity of trigonometric function, the phase value of each point in the interference image can only fall within the range of principal value (, ], i.e. wrapped phase, and the real interference phase is obtained by adding an integral multiple of 2 based on the wrapped phase, namely, the constant n in the above-mentioned calculation formula, so the real interference phase must be obtained by phase unwrapping. In the present invention, a minimum cost flow algorithm is used to perform the phase unwrapping operation. In order to ensure the quality of phase unwrapping, the periodic mean filtering method is used to filter out the interference of the noise before the phase unwrapping so as to make the interference phase clearer.

[0047] Step (4), PS point selection: during the deformation monitoring of GB-SAR system in a long period of time, as a result of the effect of temporal decorrelation, the scattering characteristics of many target points will change. Only the PS point can be used to perform appropriate APS compensation. In the present invention, the PS points are selected by using the joint amplitude and amplitude dispersion index and coherence coefficient method. First, since the PS points need to be updated, SAR images are grouped during processing, every F SAR images are set as a group, and PS candidate selection is performed for each group. Then, the correlation coefficient information of all PS candidates is used to perform the secondary screening for PS candidates, and the final set of PS points is selected. Step 4 specifically the following steps.

[0048] 4.1) PS candidate selection: the PS candidates are selected by using the amplitude and the index of the amplitude dispersion.

[0049] 4.1.1) The PS candidates are selected by using the amplitude dispersion information. The calculation formula of amplitude dispersion is as follows:

[00013]${}_{\mathrm{amp}}=\frac{{}_{\mathrm{amp}}}{m_{\mathrm{amp}}}$

[0050] where .sub.amp represents the time sequential amplitude standard deviation of the target points, and m.sub.amp represents the time sequential amplitude average of the target points. In the process of amplitude dispersion selection, the amplitude dispersion threshold .sub.adi is set first, then the pixel points having amplitude dispersion less than the threshold .sub.adi are selected as PS candidates. In the GB-SAR data processing process, generally, the amplitude dispersion threshold .sub.adi is set as 0.1.

[0051] 4.1.2) The PS candidates are further selected by using the amplitude information to further select. In SAR images, there are points having small amplitude dispersion and large phase error. Therefore, after using amplitude dispersion information to select the PS candidates, the amplitude information of the target point can be reused to select again. The amplitude selection method is expressed by a mathematical formula as follows:

[00014]$\{\begin{array}{c}{}_{\mathrm{thr}},\mathrm{PS}.Math..Math.\mathrm{candidate}\\ <{}_{\mathrm{thr}},\mathrm{reject}.Math..Math.\mathrm{point}\end{array}$

[0052] where is the average amplitude information of the P candidates selected according to the amplitude dispersion information in F images participating in the selection, and .sub.thr is the amplitude threshold. The quality of PS points can be further improved by the method of selecting again using the amplitude information.

[0053] 4.2) Secondary screening: the secondary screening is performed by using the coherence coefficient information. The coherence coefficient is an important index to evaluate the interferential phase accuracy of two SAR images. Therefore, the PS candidates are subjected to the secondary screening according to the coherence coefficient information, and the final set of PS points can be selected. The window with the size of is usually used to calculate the coherence coefficient. The mathematical expression of the coherence coefficient is as follows:

[00015]$=\frac{\underset{}{.Math.}.Math.\underset{}{.Math.}.Math.M.Math.S^{*}}{\sqrt{\underset{}{.Math.}.Math.\underset{}{.Math.}.Math.{.Math.M.Math.}^2.Math.\underset{}{.Math.}.Math.\underset{}{.Math.}.Math.{.Math.S.Math.}^2}}$

[0054] The coherence coefficients of each PS candidate in the current image and the previous image are calculated, then all points having the coherence coefficient less than the threshold .sub.coh are removed by setting the coherence coefficient threshold .sub.coh, the remaining PS candidates are the results of the final selection of PS points. In general, the coherence coefficient threshold .sub.coh of the highly-correlated PS points can be set as 0.9.

[0055] Step(5), atmospheric phase compensation: the atmospheric phase is estimated by using a two-dimensional-polynomial model based on the PS points selected in step (4). Since APS is strongly correlated spatially, a two-dimensional-polynomial model with a transformation degree of r is established:

[00016]${}_{\mathrm{atm}}\left(p,q\right)=b_{00}+\underset{t=0}{\overset{r}{.Math.}}.Math.b_{t,\left(r-t\right)}.Math.{\left(p-p_{\mathrm{ref}}\right)}^t.Math.{\left(q-q_{\mathrm{ref}}\right)}^{r-t}$

[0056] where .sub.atm is the atmospheric phase, (p,q) is the coordinate of the pixel point, (p.sub.ref,q.sub.ref) is the coordinate of the reference pixel point, and b.sub.t,(r-t) is the polynomial model parameter. The APS model parameter is estimated based on the PS points selected. Ideally, the deformation phase of the PS point selected is zero, so the following equation can be obtained:

(p,q)=.sub.atm(p,q)+.sub.noise(p,q)

[0057] In the present invention, the model parameter {circumflex over (B)}.sub.t,(r-t) is estimated by least squares method, and then the atmospheric phase of the deformation area is calculated by using the estimated model parameter {circumflex over (b)}.sub.t,(r-t). After the atmospheric phase compensation, the true phase value caused by the deformation can be obtained, so that the deformation of the monitoring area can be calculated.

EMBODIMENT

[0058] The parameters of GB-SAR system in this embodiment are set as Table 1 shown below.

TABLE-US-00001 TABLE 1 System parameters Carrier frequency 17.2 GHz Bandwidth 320 MHz Synthetic aperture length 0.85 m Synthetic aperture time 15 s Spatial sampling interval 0.005 s Range resolution 0.48 m 3 dB antenna beam width 45 degrees

[0059] The experimental site is located in Yanguo Gorge, Yongjing County, Gansu Province of China. The observation time is from 17:30 PM on Apr. 1, 2017 to 17:30 PM on Apr. 2, 2017, 24 hours for total. The deformation monitoring processing of the measured data of the mountain is completed by using the GB-InSAR-based atmospheric phase compensation method of the present invention.

[0060] FIG. 2 shows the result of the PS points selected by using amplitude dispersion information, amplitude information, coherence coefficient information, and the joint triple threshold method of the present invention, respectively. Two PS points A.sub.1, A.sub.2 are selected to analyze the time sequential evolution of the differential phase before and after the atmospheric phase compensation, as shown in FIG. 3 (a) and FIG. 3 (b), respectively. The dot curve shows the time sequential evolution of the differential phase before the atmospheric phase compensation, the diamond curve shows the result after APS compensation using a range-based linear model, and the asterisk curve shows the result after atmospheric phase compensation using the method of the present invention. As one can tell from FIG. 3, the atmospheric phase can be more accurately estimated by the method of the present invention. The conclusion can also be further verified by the unbiased estimation of the variance error shown in FIG. 4. The unbiased estimation of the variance error is defined as: S.sup.2=(.sup..sup..sub.atm)/(3), wherein represents the original interference phase vector of PS point, .sub.atm represents the estimated atmospheric phase vector of PS point, represents the number of PS points, and represents the conjugate transpose operation. As one can tell from the result of this embodiment, the method of the present invention can realize the estimation and compensation of the atmospheric phase effectively, and has certain practicability.

[0061] To conclude, the above-mentioned embodiments are merely preferred embodiments of the present invention and not intended to limit the present invention. Any modification, equivalent substitution, and improvement derived without departing from the spirit and principle of the present invention shall be considered to fall within the scope of the present invention.