Three-Dimensional Ground-Based SAR Imaging Algorithm Based on Keystone Formatting and Subblock

Three-dimensional (3D) ground-based synthetic aperture radar (GB-SAR) draws attention because of the ability to obtain high accuracy 3D information of the monitoring terrain. Besides, the 3D GB-SAR system is welcome due to the flexibility and large coverage in angle and range. However, existing 3D imaging algorithms for GB-SAR data focusing have limits of high computational complexity or narrow applicable scope. To realize 3D displacement monitoring with high spatial resolution and short revisit time, in this paper, a novel 3D imaging algorithm is proposed. Based on characteristics of the model of echo data from the large range and wide-view angle scenario, the proposed method uses keystone formatting to complete range migration correction and subblocks dechirping to realize horizontal focus. With the method, the reflectivity of the monitoring scenario is obtained by only one time of linear interpolation and several times of fast Fourier transforms. The main advantages of this algorithm are high imaging precision and low computational cost, and in addition, it is applicable for large illuminating coverage including the near-field and the far-field of the radar aperture. The imaging results are sampled on pseudo-spherical grid, aiming to simplify the formulation. Finally, this method is extensively validated with simulations.


I. INTRODUCTION
Ground-based synthetic aperture radar (GB-SAR), a newly rising synthetic aperture radar technology, can realize large-area, all-weather and all-time monitoring with high accuracy and short revisit time. In addition, the GB-SAR system applies the differential interference (DIn) techniques of space-borne SAR to catch up the deformation in millimeter or even sub-millimeter precision. Therefore, it is widely used for deformation monitoring of open pit mines [1], [2], landslides [3], [4], buildings and structures [5] and glacier snow mountains [6], [7].
The conventional two-dimensional (2D) GB-SAR imaging system can only acquire an image which actually is The associate editor coordinating the review of this manuscript and approving it for publication was Gerardo Di Martino . a two-dimensional projection on the range-azimuth plane of the natural three-dimensional space. Therefore, the 3D GB-SAR system is adopted to obtain targets' height information. It obtains two-dimensional large aperture through the radar motion in azimuth and height directions which enables it to have azimuth and height resolutions, and further realizes 3D high precision imaging combined with transmission of the large time-bandwidth signal. In practical application, DInSAR system is used to obtain the tiny deformation of the 3D monitored terrain. The key step of the processing is the registration of the radar image and 3D DEM (Digital Elevation Model) [8]. That is corresponding the pixel points of the radar image to the points in the 3D DEM. And the deformation information of the 3D DEM can be obtained by the deformation of the corresponding pixel points. However, the 2D radar image can only obtain the 2D information of the 3D monitored terrain, so only the 2D deformation information of 3D terrain can be gotten. On the contrast, after the registration of the 3D radar image and 3D DEM, every point of the 3D DEM can be match to the pixel points of the 3D radar image. Based on the 3D GB-DInSAR technique, the 3D deformation information of the pixel points can be obtained, so the 3D displacement information of 3D terrain can be gotten without the information loss. In addition, 3D GB-DInSAR system can be used to achieve the 3D deformation monitoring of snow mountains, landslides, volcanoes or manmade structures. The tomographic 3D GB-SAR system can be used to measure the high-dimensional information of snow mountains. For example, the SAPHIR team utilized a tomographic 3D GB-SAR system transmitting signals in X-band and Ku-band to obtain the vertical structure of snow mountains and such natural environment [9]. Stephen Joseph Preston and David G. Long made use of 3D GB-SAR to monitor snowpacks [10], [11].
The key issue of 3D GB-SAR system is the fast and accurate imaging algorithm. The echo data produce the 3D radar image after imaging processing, so the processing computational cost directly determines the system's real-time performance. The back projection algorithm (BPA) is a widely used algorithm because of its high precision and simplicity. In addition, it can be easily used in SAR systems with all kinds of configurations. However, BPA's high computational cost makes it not suitable for large-scale and fast-deforming monitoring. Even though there are some proposed methods to improve BPA's performance [12], [13], it still cannot meet the requirements of dealing with the fast-changing phenomenon. Reference [14] has proposed a far-field pseudopolar format algorithm (FPFA). It only retains one-term in the range history of targets. Therefore, FPFA is only valid in near-field even if it is fast to focus data. There is a wavenumber domain imaging algorithm, RMA [15] It obtains the spatial spectrum of targets from the echo data, then combined with deformation, filling and etc. of the spatial spectrum to obtain the scattering function of the whole scene. However, the stolt interpolation process of RMA is a necessary step but has high computational load [16], [17]. Thus, RMA is not applicable for real-time imaging of 3D GB-SAR. This paper proposes an efficient and accurate 3D imaging algorithm, that is Three Dimensional Keystone formatting and Subblock Dechirping (3D-KSD) algorithm. The core of this algorithm is the range cell migration correction (RCMC) and horizontal dechirping processing [18], [19], where the RCMC is realized by keystone formatting, and the horizontal focusing can be achieved by subblocks dechirping operation. The approximate preserving quadratic terms of the range history makes 3D-KSD algorithm suitable for near-field and far-field focusing. At the same time, because of the only once keystone formatting and several fast Fourier transforms (FFT), it has a low computational cost. Besides, the number of horizontal subblocks decreases with the range increasing, which reduces the computational cost further. Thus, the algorithm can achieve fast focus of a wide and large scene. This paper is organized as follows. In Section II, the principle of three-dimensional imaging radar is elaborated and the mathematical model of pulse-compressed echo signal is given, based on which the imaging algorithm is presented in detail. In addition, the application conditions and computational complexity of the algorithm are discussed. In Section III, the validation of the algorithm is extensively summarized with the simulated data and a comparison with BPA and FPFA is illustrated. Finally, our conclusions are outlined in Sections IV.

II. THE BASIS OF THE THREE-DIMENSIONAL GROUND-BASED SAR SYSTEM
The 3D GB-SAR system achieves cross-range resolutions in azimuth and height directions by a 2D synthetic aperture. Here, the fuzzy function is applied to analyze the resolution performance of the 3D GB-SAR system, whose data collection geometry is shown in FIGURE 1. The antenna moves in the track as the arrow along the x direction shows. When the antenna arrives at the end of the track, the track moves up the distance between two sample points as the arrow along the y direction shows. Thus, we can see that the length of the track is L which is the synthetic aperture in the azimuth direction, and the distance that the track moves is H which is the synthetic aperture in the height direction. The movement of antenna in x−y plane achieves two-dimensional synthetic aperture. The colored dots on the x−y plane are the sampled points of twodimensional synthetic aperture where the antenna transmits and receives signal and correspondingly the antenna location is P(n, i) = (x n , y i , 0). The three-dimensional area enclosed by the purple dotted line is the detection area of the 3D GB-SAR system. There is one point target P 0 (ρ 0 , θ 0 , φ 0 ) in the area. To avoid data redundancy, the imaging result is formed in a pseudo-spherical coordinate where the target's location is P 0 (ρ 0 , θ 0 , φ 0 ). Different with spherical coordinate, the azimuth angle θ 0 is the angle between the OP 0 vector and the y-z plane, and the height angle φ 0 is the one between the vector and the x-z plane. ρ 0 is the distance from P 0 to the aperture center O (|OP 0 |) shown as the right picture in FIGURE 1. VOLUME 8, 2020 Transmitting signal of the system is chirp signal, which is s t (t) = rect t T p exp(j2πf c t)·exp(jπKt 2 ). t is the fast time; T p is the pulse width; f c is the carrier frequency; K is the chirp rate. And the range history from the target P 0 (ρ 0 , θ 0 , φ 0 ) to the antenna P(n, i) is: R(n, i; P 0 ) = x 2 n +y 2 i +ρ 2 0 −2x n ·ρ 0 sin θ 0 −2y i ·ρ 0 sin φ 0 (1) Thus, the echo signal from P 0 is: where τ n,i = 2·R(n, i; P 0 ) c is the time delay of target P 0 . The fuzzy function was first proposed by J.Ville, and has become an effective mathematical tool for studying radar signal and getting its resolution performance. Here, we use the fuzzy function to analyze the range resolution, azimuth resolution and height resolution of 3D GB-SAR system. The fuzzy function is defined by the correlation operation between the echo signal and the transmitting signal [20], [21]. So it is written as, (3) shown at the bottom of the next page.
where the integral term is the range fuzzy function χ R (τ n,i ). The range resolution is determined by where A r is the delay resolution constant; U (f ) is amplitude spectrum of the transmitting signal. And (3) can be written as: Since the point P 0 is in the far field of the system, the τ n,i is as follows using Taylor expansion at P 0 : Thus, (5) can now be expressed as: where χ A ( R n,i ) is the azimuth fuzzy function, and χ H ( R n,i ) is the height fuzzy function of the system. What can be concluded is the system's fuzzy function can be divided into three fuzzy functions along the range, azimuth and height directions under the far field. So, according to (4)and (7), the 3D resolutions respectively are: where B is the bandwidth of the transmitting signal; 0.886 is the width factor of the -3dB beam width; ρ is the target range; δ r is the range resolution; δ a and δ h are the azimuth and height resolutions.
The main characteristics of the 3D GB-SAR system as shown in FIGURE 1 are summarized as follows: 1) Short aperture lengths in azimuth and height directions (1-4m): angle resolutions (-3dB beam width) are constant and space resolutions in cross-range directions (δ a and δ h ) are increasing linearly with the increase of the range, as shown in (8). 2) Large antenna beam width in azimuth and height directions in order to illuminate the whole monitoring scene. 3) A large range scope: targets in near-field and far-field of the system are covered by the antenna illuminating. The baseband range-compressed echo signal from the point P 0 is where p r (t) is the compressed range envelope. The range history R(n, i; P 0 ) is a key point affecting the range migration and azimuth and height modulation, which can be described as follows using the Taylor expansion at the aperture center (x n = 0, y i = 0): Here, the higher order terms are contained in o r4 (·) which is less significant given L, H ρ. To further simplify the echo signal, we analyze the range history R(n, i; P 0 ) of the range envelop term and the exponential order term in (9), to decide which term of R(n, i; P 0 ) can be omitted under the system parameters in TABLE 1.
First, according to (10), the range cell migration R in the range envelop term is, where R linear is linear range migration, and R curvature is range curvature. If their absolute values are less than a quarter of the range resolution, they can be ignored. According to the system parameters, the maximum of R linear is, Under the system parameters in TABLE 1, it cannot always less than a quarter of the range resolution and should be considered.
The maximum of R curvature is, That is Next, the phase of echo signal in (9) is When the absolute value of phase term is less than π 8 , it can be ignored. The second-order term of the phase is ψ quadratic n, i; P 0 When we limit ψ quadratic n, i; P 0 max = πL 2 λ c ρ min ≤ π 8 , we can get ρ min ≥ 8· L 2 λ c . It cannot be satisfied all the time in 3D GB-SAR system. So, the second-order ψ quadratic n, i; P 0 cannot be ignored.
The third-order term of the phase is ψ cubic n, i; P 0 Therefore, when the range ρ satisfies (15) and (19), we can retain R(n, i; P 0 ) to the first-order of the range envelop term and the second-order of the exponential order term of the echo signal model in (9), ignoring the range curvature term in the range envelop term and high order terms in the phase. Consequently, it can be rewritten as: which is the base for following algorithm formulation. The simplified signal model demonstrates that the range migration of the echo signal is linear and the phase history is parabolic both in azimuth and height directions. Thus, the 3D-KSD algorithm is proposed. FIGURE 2 is the flowchart of the algorithm.

III. RANGE CELL MIGRATION CORRECTION
In order to focus the energy of the point P 0 on its range ρ 0 , we apply Keystone formatting, which is a useful tool in SAR moving target processing [22], to realize RCMC. Firstly, (20)after Fourier transformation in range is where f is range baseband frequency; P r (·) donates range envelope in frequency domain. Next, replace the azimuth variable x n and the height variable y i with f c f +f c x m and f c f +f c y j for each range frequency, so (21) can be rewritten as 90122 VOLUME 8, 2020 where x m and y j are respectively the new azimuth variable and the new height variable after Keystone formatting. Lastly, the free of migration time signal after inverse-range Fourier transformation is described as Since B r /f c 1 in 3D GB-SAR systems, we can simplify rect (23)

IV. AZIMUTH AND HEIGHT FOCUSING
The second part in FIGURE 2 aims to accomplish the azimuth and the height focusing. Firstly, equation (23) should be cross-range transformed into range-Doppler domain, that is, pseudo-spherical domain. In addition, the x m 's Fourier counterpart has a linear relation with sin θ which is f m = 2 sin θ/λ c [22]; the y j 's Fourier counterpart has a linear relation with sin φ which is f j = 2 sin φ/λ c . Therefore, the transformed signal can be described in the ρ−sin θ−sin φ domain. According to the principle of stationary phase, the relation between x m and sin θ can be written as Similarly, the relation between y j and sin φ is where K a and K p , the frequency modulation slopes of the azimuth and the height, are as following, respectively.
Thus, the transformed signal in the ρ−sin θ−sin φ domain is where L sin θ (ρ 0 , θ 0 ) is the width of support field in sin θ domain, and L sin φ (ρ 0 , φ 0 ) is the width of support field in sin φ domain We can observe the equation (23) which shows that different targets have the same azimuth support from −L/2 to L/2 and the same height support from −H /2 to H /2 in ρ−x m −y j domain. Whereas, equation (27) shows that targets have segregated supports in the ρ−sin θ−sin φ domain, centering on their angle coordinates (sin θ 0 , sin φ 0 ) with the azimuth extension L sin θ (ρ 0 , θ 0 ) and the height extension L sin φ (ρ 0 , ϕ 0 ).  FIGURE 5 (c). We can see that energy of different point targets after cross-range transformation is dramatically separated. Compared with the length of sin θ and sin φ axis, the targets' energy support is small. The blocking strategy as following is based on the energy-separation character of the 3D GB-SAR signal. VOLUME 8, 2020

A. AZIMUTH AND HEIGHT BLOCKING
Besides Keystone formatting, another key point in 3D-KSD algorithm is the cross-range (azimuth and height directions) blocking, by which targets in monitoring scene are segmented into several subblocks. The purpose of blocking is to realize azimuth and height focusing. Because of blocking, spatially variant dechirping can be applied to focus the targets with different K a and K p . The principle of blocking is given below at first and the derivation would be elaborated in detail later. Given a range gate ρ, the maximum width of blocks is where |sin θ+sin φ| max is the maximum absolute value of the sum of sin θ and sin φ, and sin 2 θ+sin 2 φ max is the maximum of the sum of sin 2 θ and sin 2 φ. In order to simplify the blocking process, we limit the blocking width in azimuth direction is equal to that in height direction, that is, sin θ (ρ) = sin ϕ (ρ) = (ρ). An example of 3D GB-SAR system parameters is given in TABLE 1, which is used for following derivation.  Under the parameters in TABLE 1, FIGURE 6 (a) illustrates the variation of (ρ) with the variable ρ which shows that (ρ) increases as range ρ becomes larger. The minimum value of (ρ) is 0.0287 at the nearest range of 100m. When (ρ) is larger than the length of sin θ or sin φ axis, this dimension would not be divided into blocks. Considering effectiveness of computation, the width of blocks should be equal to the upper limit of (ρ) in (29). Now, the blocking processing in azimuth direction is shown as following and illustrated in FIGURE 7.
1) At the range gate ρ, the number of sample points of each normal subblock is where υ sin θ is the sampling interval in sin θ axis, and the concept of the normal subblocks means blocks with the same number of points. N a is the number of sample points of sin θ axis, and is to round down numbers.
2) The number of the subblocks in range gate ρ is, where is round up to numbers. 3) If N a can be divided by n (ρ) exactly, the N a sample points of the range ρ can be uniformed segmented into m (ρ) subblocks each with the number of points n (ρ), as shown in the upper figure of FIGURE 7(a). If N a can be divided by n (ρ) exactly with a remainder, the N a sample points of the range ρ still will be segmented into m (ρ) subblocks. The m (ρ)−2 subblocks have the same sample points n (ρ), located in the middle of the sin θ axis. The other two subblocks are located in the end of the sin θ axis respectively with sample points n 1 (ρ) and n 2 (ρ), as shown in the lower figure of FIGURE 7 (a). n 1 (ρ) and n 2 (ρ) are two near integers and satisfy the next equation that is n 1 (ρ)+ n 2 (ρ) = n (ρ)+rem (N a /n (ρ)).
Following the instructions above are the same with the azimuth blocking processing in the reference [22] which presents the 2D-KSD algorithm. The results of azimuth blocking is shown in FIGURE 6 (b) and (c) according to the TABLE 1. In the simulation, the number of azimuth sample points is N a = 512, and the sample interval in sin θ axis is υ sin θ = 0.0045. Firstly, the length of the normal subblocks computed by (30) increases as the range gate ρ is bigger as illustrated in FIGURE 6(b). Then the number of subblocks m (ρ) computed by (31) is shown in FIGURE 6(c). To ensure high accuracy, the nearer the range ρ is, the more the subblocks m (ρ) are. At the nearest range of 100m, m (ρ) is 53, and with the range ρ increasing, m (ρ) decreases until it drops to 1 when the length of a subblock n (ρ) is equal to the length of sin θ axis. As a consequence, tighter blocking in near range and looser blocking in far range indicate the computational complexity descending, as shown in FIGURE 7(b). After m (ρ) and n (ρ) are obtained, step 3) and 4) give the procedures of the data segmentation of each range ρ. The height blocking instructions are the same with the azimuth and the number of subblocks is p (ρ).
The blocking process can be achieved mathematically, that is, the (i, k) subblock can be obtained by the product of data and window function W (sin θ, sin φ; ρ, i, k) ,i.e., where i and k mean the subblock is the ith in azimuth and the kth in height. The subscript '0' indicates the parameters about targets. The horizontal window function to obtain the horizontal subblock is which means that the subblock's azimuth center, sin θ i , azimuth length, L i sin θ , height center, sin φ k and height length, L k sin φ , all vary with the range ρ. According to (32), the window function (33) segments the data support domain rect sin θ −sin θ 0 , as shown in FIGURE 7 (c). Reversely, connect all the segments and gain original data support profile.
Transform the m (ρ)·p(ρ) subblocks to x m −y j domain and get Connect all segments and obtain the original data profile, i.e., Similarly, the segmentation in sin φ domain means the segmentation in y j domain, and the y k 0 and L k y can be written as All segments in y j domain are connected to gain original data profile, i.e., Then all subblocks are dechirped by multiplying with reference signal in ρ−x m −y j domain. The different reference signal is determined by the (i, k) subblock's center (sin θ i , sin φ k ) and the range ρ as below.
Adding all the product of (35) and (42) at range gate ρ, the dechirped result is where the second line is what we want. The last line is residual phase error, denoted by ε de ρ, x m , y j ; ρ 0 , θ 0 , φ 0 , introduced by the difference between the angle (θ i , φ k ) of the reference signal and the angle (θ 0 , φ 0 ) of the target. Large value of ε de ρ, x m , y j ; ρ 0 , θ 0 , φ 0 would result in defocus or low phase accuracy. The denser the blocking is, which means (sin θ i , sin φ k ) is closer to the real value (sin θ 0 , sin φ 0 ), the smaller ε de ρ, x m , y j ; ρ 0 , θ 0 , φ 0 is. But denser blocking brings in large computation cost. Thus, in order to balance the imaging quality and the computational complexity, the 3D-KSD algorithm will take measures to control the phase error. Firstly, we limit the value of (ρ) mentioned above to ensure ε de ρ, x m , y j ; ρ 0 , θ 0 , φ 0 ≤ π/8. Here, the phase error in range ρ 0 is According to (38) and (41), each subblock has different and mutually exclusive data profiles, and the connection of the data profiles has consistent amplitude rect x m L rect y j H . So we only need to limit the phase 2π λ c ρ sin θ i x m +sin φ k y j 2 − sin θ 0 x m +sin φ 0 y j 2 to less than π/8 for any i and k. 2π where L sin θ and L sin φ are the data support in the azimuth and height direction. The 3D GB-SAR system commonly has the same azimuth resolution and height resolution, so the lengths of the synthetic aperture in azimuth and height are the same (L = H ). , the width of subblocks, needs to be determined ( sin θ = sin φ = ). Under the limit of π/8, the maximum of is According the parameters in TABLE 1, the blocking result up of the whole scene is shown in FIGURE 8. It illustrates that up increases with the increase of range ρ 0 and the decrease of |sin θ 0 +sin φ 0 |. In addition, the variation of range has more influence than the angle on up . Here, in 3D-KSD, the minimum value up is chosen as the blocking width (ρ) in each range, i.e., (ρ) can satisfy the block width limit of all targets with different locations (θ 0 , φ 0 ) in range ρ. It is exactly the blocking width given in (29). Under the limit of blocking width, the horizontal focus can be achieved by transforming (43) into sin θ−sin φ domain.
where the profile is the ideal focus result. The limit of phase error π/8 is enough to meet the high phase precision requirement on most occasions.

C. APPLICABLE SCOPE AND COMPUTATIONAL COMPLEXITY
In this section, the comparison between 3D-KSD and other three 3D imaging algorithms (BPA and FPFA) is given below. BPA has high imaging accuracy and can focus data both in the far-field and near-field, but high computational cost makes it not suitable for monitoring fast deforming targets. RMA also can obtain high-quality imaging results, but has strict requirements for the radar system's configuration as well as high computational cost. FPFA can only focus the far-field targets efficiently.
3D-KSD's applicable scope is determined by (15) and (19) i.e., where max {} indicates the maximum value of the input variable. We can see that the minimum range ρ min and the angle scope |sin θ+sin φ| max of the applicable scope are mutually restricted, as shown in FIGURE 9, and four sets of typical scene parameters are given in TABLE 2 according to (49). The larger the range is, the narrower the horizontal scene is.  On the contrast, FPFA has a limited applicable scope given by indicating the minimum range of imaging scope is 864m according to TABLE 1. So FPFA cannot obtain the near field imaging. Then, the comparison on computational complexity of 3D-KSD and BPA is given using floating point operations (FLOPs) [18]. Assuming that the original dataset (after range compressed) is N r ×N a ×N p , where N r is the number of range sample points, N a is the number of azimuth sample points, and N p is the number of height sample points. Given that the interpolation kernel length is M ker , and the number of subblocks in azimuth is m (ρ), in height is p (ρ). The calculation results of 3D-KSD's FLOPs is given in TABLE 3. Assuming that in the dataset N a = N p = N r = N , the computational complexity of 3D-KSD is O N 3 log (N ) . However, the computational complexity of BPA is O N 4 . Moreover, some existing methods based on FFT and complex multiplications are used to achieve Keystone formatting to reduce the computation further. It should be noted that, when the range increases, the number of subblocks decreases, resulting in the dramatically reduce of 3D-KSD's computation. In far field, the linear terms in R(n, i; P 0 ) (10) is enough for focusing, so the 3D-KSD becomes to FPFA. At this time, the horizontal blocking is not needed, and the reference signal is 1. The focus in horizon can be achieved by FFT after RCMC, which makes 3D-KSD more efficient than BPA.

V. SIMULATIONS
To verify the feasibility of 3D-KSD and obtain its advantages compared with other 3D imaging algorithms, the imaging simulation experiments are given in this section. The system parameters are shown in TABLE 1.
First, the imaging process of the point target located in (60m, 30 • , 30 • ) is given in FIGURE 10. FIGURE 10(a) gives the range-compressed signal before Keystone formatting, showing the obvious range cell migration in range-azimuth domain. FIGURE 10(b) shows that after Keystone formatting, the range cell migration is corrected and the target energy is totally centralized in the target range, 60 m. Then the data is transformed into ρ−sin θ−sin φ domain. After data segmentation, the subblocks are obtained and inversely transformed into ρ−x m −y i domain. By multiplying subblocks with reference signals and summing subblocks in each range, the dechirping operation is finished. After the data being inversely Fourier transformed into ρ−sin θ−sin φ domain, the horizontal focus is achieved as shown in FIGURE 10(c) and (d). The 3D imaging slice of the point target and onedimensional imaging slice results along the range, azimuth and height are shown in FIGURE 11. We can see that the target is well focused and the performances of the three directions are dramatically good. We can see from FIGURE 13 (b), the side view and (d), the front view of the 3D imaging result that 3D-KSD has excellent focusing performance. In addition, the positions of focusing points in 3D imaging result are totally consistent with the targets position we set. This experiment validates that 3D-KSD is competent for multi-target imaging. Thus, 3D-KSD is qualified for imaging complex terrain.
Next, the comparison between 3D-KSD, BPA and FPFA is given by simulation data. First, there is a target at (60m, 0 • , 0 • ), locating in the near-field of the radar aperture.    three algorithms. We can conclude that 3D-KSD has the same performance with BPA, and they are both valid in near-field. However, FPFA cannot focus the target. Then, FIGURE 14 (d), (e), and (f) show the imaging results of the target at (500m, 0 • , 0 • ), in the far-field of the radar aperture. The target is well focused by the three algorithms. We can conclude that the three imaging algorithms have good performance in the case of far-field and the quantitative values are given in the TABLE 4.
In addition, the computation time required by imaging algorithms is given in TABLE 5. The time cost comparison of 3D-KSD, BPA and FPFA demonstrate the 3D-KSD has high computation efficiency in near-field and far-field while ensuring good focusing performance.
Finally, the performance comparison of 3D-KSD, BPA and FPFA is concluded. An imaging algorithm's performance includes focusing accuracy and computational complexity.  So VI gives the comparison result of imaging algorithms based on the two aspects in near-field and far-field. VI shows that, in far-field and near-field, 3D-KSD has lower computational complexity and the same high focusing accuracy with BPA. In far-field, FPFA is good in the two aspects. However, in near-field, FPFA cannot focus, so it cannot be used in the situation where both near-field and far-field imaging are required. In conclusion, the proposed method, 3D-KSD, is suitable for large range scope, high precision and real time imaging.

VI. CONCLUSION
This paper proposes a 3D imaging algorithm using for fast and accurate imaging of 3D GB-SAR system. It can be widely used for 3D GB-SAR imaging systems and most of 3D radar imaging systems with similar configurations and radar parameters. This method driven from the analysis and modeling of the echo data of 3D GB-SAR system is an algorithm based on Keystone formatting and Subblock Dechirping. In terms of RCMC, the Keystone formatting is used to eliminate the linear component of the echo data. Then the algorithm realizes the horizontal focus by subblock dechirping. The whole imaging process includes linear interpolation, FFTs and complex multiplications which dramatically save the computation cost. In addition, the 3D-KSD algorithm can be used in both near-field and far-field effectively. Thus, 3D-KSD can be used to monitor large range scope filed and obtain its 3D information. In addition, low computational complexity makes the proposed imaging algorithm fast enough to obtain the 3D tiny deformation of speedy deforming targets such as landslides, buildings and structures and glacier snow mountains. The 3D GB-SAR system equipped with 3D-KSD is a powerful tool to predict geological disasters caused by tiny deformation.