Tomographic Imaging for Orbital Radar Sounding of Earth's Ice Sheets

BingSat-Tomographic Observation of Polar Ice Sheets (BingSat-TOPIS) is a spaceborne multistatic radar sounding system that can achieve high resolution and stereoscopic observation. It is designed to penetrate ice sheets and acquire tomographic image. The satellite group fly over the Polar Regions, which can form about 6.56-km cross-track baseline. This cross-track baseline is formed by one master satellite and 40 slave CubeSats. In this article, we propose a tomographic imaging algorithm for orbital radar sounding of the Earth's ice sheets. First, we give the method to calculate the two-way slant range in air and ice. The phase errors of the method are smaller than π/4 rad. Second, we express the radar data cube for a point target in ice sheets. The radar data cube is made of 40 bistatic synthetic aperture radar echo signals. At last, the echoes are simulated for the TOPIS system, and a matched filter is used for pulse compression in range. Then, the back projection algorithm is used in track and cross-track imaging. The experimental results demonstrate that our tomographic imaging algorithm is effective and reliable.


I. INTRODUCTION
O NE of the major purposes of the polar observation is to gather the information of ice mass in Polar Regions, which has an important influence on global climate change and sea level rise [1], [2], [3]. There are many kinds of radar sounding systems to monitor the losing mass of polar ice sheets [4], [5], [6], [7], [8]. Currently, vehicle-borne and airborne systems are the most commonly used platforms. They cost a lot of manpower and material resources. Besides, the update rate in time is still not enough. In the northern part of Greenland, it is reported that the sum of attenuation is more than 60 dB for 3-km thick interior ice [9]. Therefore, the radar sounding system called BingSat-Tomographic Observation of Polar Ice Sheets (BingSat-TOPISs) is proposed to achieve high resolution and stereoscopic observation [10]. It can get fine resolution in cross-track, along-track (track or azimuth), and vertical (range). Manuscript  For BingSat-TOPIS, it faces a lot of challenges [11]. This article is mainly focused on the tomographic imaging algorithm [12] for the internal layer in ice sheets, which is not just simple combination of 2-D (track × vertical) images dataset. BingSat-TOPIS has the following characteristics.
2) The received signals travel through two media (air and ice). In this article, we mainly propose a tomographic imaging algorithm for orbital radar sounding of Earth's ice sheets. In Section II, system formation of BingSat-TOPIS and some key parameters are described. Based on Earth observation geometry, the slant range and the signal model for BingSat-TOPIS are built in Section III. In this part, we give the equivalent slant range of a new model to improve computational efficiency, which only needs to calculate the quintic equation twice for a two-way slant range. The new model is used to reduce the computation time. The detailed description of simulation processing and tomographic imaging results is presented in Section IV. Finally, Section V concludes this article.

II. BRIEF INTRODUCTION FOR THE BINGSAT-TOPIS SYSTEM
BingSat-TOPIS is designed to realize a 100 m × 20 m × 5 m resolution (cross-track × track × vertical), which can compensate 65-dB ice attenuation and realize -3.5-dB system sensitivity [10]. It has sounded the ice depth about 4000 m. The system can realize the high-resolution tomographic imaging of polar ice sheets by combining various advanced technologies.
The vertical resolution ρ v for the nadir-looking radar sounding system is given by where c ice is the speed of the electromagnetic wave in ice, c is the speed of light in vacuum, and B r is the bandwidth of the transmitting signal. The ice surface resolution ρ a is defined as where D a is the length of the azimuth antenna. The length of a 40-m diameter parabolic transmitting antenna is 40 m.
In the wavenumber domain, the cross-track spatial resolution ρ c can be expressed as where φ i is the range of the incident angle, and φ e is the range of the emergence angle. In the single-input and multiple-output (SIMO) formation, the range of the incident angle φ i is zero, and the range of the emergence angle φ e depends on the maximum cross-track baseline B ⊥ . Therefore, ρ c becomes Table I gives the system parameters of BingSat-TOPIS [10]. As shown in Fig. 1(a), the master satellite in the center of the formation carries a large transmitting antenna, and the 40 slave satellites receive the echoes together and retransmit them to the master satellite for sampling and recording. In Fig. 1(b), in order to reach the Polar Regions simultaneously and keep a stable baseline, all satellites must fly in nearly circular with only small differences in eccentricity.

A. Slant Range Between Air and Ice
BingSat-TOPIS is different with the ground penetrating radar (GPR) system and airborne radar sounding system. For the satellite case, the rectilinear geometry is replaced by the curved Earth geometry. The orbits of the master and the slave satellites are curved. The ice sheets on Earth are also curved. Besides, the Earth is rotating independently of the satellite orbit.
For BingSat-TOPIS, the master satellite in the middle is surrounded by 40 slave satellites. It carries a high transmitting power antenna, which is vertical observation. The slave satellites each carries a dipole receiving antenna and receives the radar echo from Earth. In this article, unless otherwise stated, all parameters are expressed in Earth centered rotating (ECR) coordinates.
BingSat-TOPIS is composed of 40 bistatic SAR configurations. A bistatic SAR configuration is chosen and analyzed. First, an accurate equation of slant range is deduced. Based on the accurate equation, a quintic equation for one-way slant range is presented with shorter computation time than the accurate equation.
Geometry of slant range (one-way) between antenna and a point target in ice sheets is illustrated in Fig. 2. As shown in Fig. 2, S represents the position of the satellite, ϕ is the crosstrack beamwidth, h is the height of the satellite, O is set as the Earth's center, and R L is the local radius of the Earth. There is a target P 2 , whose ice depth is d. P 1 is the incident point at air-ice interface between the satellite S and the target P 2 . θ i1 and θ t1 are the angle of incidence and the angle of refraction at P 1 , respectively. θ i2 is the angle of incidence at P 2 . α 1 and α 2 are the angles for ∠SOP 1 and ∠P 1 OP 2 , respectively. ε 1 and ε 2 are the dielectric constant of air and the dielectric constant of ice, respectively.
Assuming t i is the azimuth time, R s (t i ) = [R sxi , R syi , R szi ] and P 2 = [P 2x , P 2y , P 2z ] are the position vectors for S and P 2 , respectively. The local radius R L is given by [13] where R e and R p are the semimajor axis and the semiminor axis of the Earth, respectively. The geocentric latitude Φ s is given by The distance between the satellite S and the Earth's center O is given by The distance between the target P 2 and the Earth's center O is given by Considering microwave propagates in air, the slant range between the satellite S and the target P 2 is given by (9) α is the angle for ∠SOP 2 . According to the cosine law, α is given by Considering microwave propagates between air and ice, the slant range is divided into slant range of air and slant range of ice. As shown in Fig. 2, the angle of incidence θ i1 is expressed as where the denominator is the distance between S and P 1 . Similarly, the angle of incidence θ i2 can be expressed as where the denominator is the distance between P 1 and P 2 . The relationship between the angle of incidence θ i2 and the angle of refraction θ t1 is given by Substituting (12) into (13), sin θ t1 is written as According to Snell's law, the refractive index of the ice layer can be written as Here, the aforementioned equation (15) is an accurate equation of slant range between air and ice. α 2 is much smaller than α 1 . If α 2 is ignored, (15) can be rewritten as where Squaring both sides of (16), we obtain If the ice depth d = 0, a quintic equation is expressed as where the solution x can be obtained by solving the quintic equation. Thus, α 2 is easily obtained through x.
For BingSat-TOPIS, it involves only a very small error between the quintic (18) shown at the bottom of this page and the accurate (15). However, the computation time using (18) is much faster than using (15).
However, the quintic (18) and the accurate (15) are both very time consuming. We need to find a way to solve this problem.
The slant range errors for a bistatic SAR in BingSat-TOPIS are shown in Fig. 3. The receiving satellite of the bistatic SAR has the farthest distance to the master satellite. As shown in Fig. 3, the maximum slant range error of the bistatic SAR in BingSat-TOPIS is no more than 0.001 m.

B. Equivalent Slant Range by Rectilinear Geometry
As shown in Fig. 4, the blue dots Sand S c represent the position of the satellite at the azimuth time t i and t c = 0, respectively. The red dot P represents the position of the point target. S c , S surf , S d , S v , P c , P surf , P , and P v are on the zero Doppler plane. The lines of P c S c , P surf S surf , P S d , and P v S v are perpendicular to the line of S c S v , and their perpendicular distances are the same. h is the height of satellite, d is the ice depth, d eq is the equivalent virtual depth, and v is called effective radar velocity. t c = 0 is the time of closest approach. t i is the azimuth time referenced to the time of closest approach. At the azimuth time t c = 0, the equivalent slant range is the sum of slant range in air R air (0, d) and slant range in ice R ice (0, d). At the azimuth time t c = 0 , l c is distance between the point of incidence and the point of refraction on the ice surface, and z c is distance between the point of incidence and the nadir of the satellite on the ice surface. Similarly, we define , l, and z. At the azimuth time t i , θ i is the angle of incidence, and θ t is the angle of refraction. At the azimuth time t c = 0, Z c = l c + z c . At the azimuth time t i , the slant range between a satellite and a point target in ice is given by As shown in Fig. 4, R air (t i , d) and R ice (t i , d) are given by where x and l vary with the azimuth time t i . The first derivative of the slant range is given by If Z c = 0, we obtain [14] dR (23) The second derivative of the slant range is written as .
The second derivative of the slant range is written as For the rectilinear geometry, the height of satellite h and the ice depth d are both constants. The beam footprint on Earth's surface for the master satellite can be written as Substituting λ = 1m, h = 449km, and L = 40m into (27), we obtain where Z c,max is about 5 km for the master satellite. The master satellite is vertical observation (449 km), the ice depth is 4 km, and the maximum cross-track Z c,max is about 5 km. There are two ways to calculate the equivalent slant range by using the rectilinear geometry of Fig. 4.
1) The equivalent slant range between a satellite and a virtual point P v in air where ) is the closest approach. The second derivative of the equivalent slant range is written as Comparing (26)  v eq t i l eq , and l eq is varies with v eq t i . The equivalent slant range in air is written as where v eq , l eq , and φ squ are the equivalent velocity, equivalent length, and equivalent squint angle, respectively.
The equivalent slant range in ice is written as According to (9), v eq is given by where t i_start and t i_stop are the start and stop time of azimuth. The equivalent length l eq is given by · v eq t i . (34) The equivalent squint angle φ squ is given by In Fig. 4, the two-way equivalent slant range errors between the bistatic SAR and the five targets in ice are calculated in two ways. The bistatic SAR is composed by the master satellite and the furthest slave satellite. Referenced to the line between the Earth's center and the master satellite, the cross-track positions of five targets are 5000 m.
As shown in Fig. 5(a), five virtual points in air and five real points in ice have the same closest approaches. Comparing with the quintic equation, their two-way slant range errors at five depths (400 m, 1200m, 2000 m, 2800 m, and 3600 m) are shown. At 2800-m depth, its two-way slant range error reaches 0.63 m, which is smaller than λ/16(0.0625 m) [15]. At 3600-m depth, its

C. Bistatic SAR Signal and Radar Data Cube
As shown in Fig. 6, the bistatic SAR geometry is given. Comparing to Fig. 5, Fig. 6 adds a receiving satellite Q.
Considering that the ice is solid ice, which has a homogeneous refractive index. At the azimuth time t i , the receiving satellite named nk, the two-way slant range of the bistatic SAR can be expressed as Rair (t i , d), and R nk−Rice (t i , d) are the distances of line SP 1 , P 1 P 2 , QQ 1 , and Q 1 P 2 , respectively. The demodulated baseband signal of the bistatic SAR (the receiving satellite named nk) for a point target in ice can be written as where the coefficient A 1 is a complex constant; w r (·) and w a (·) are envelopes of range and azimuth, respectively; τ and t i are times of range and azimuth, respectively; t c is the beam center offset time;f 0 is the radar center frequency; c is the speed of light in free space; K r is the frequency-modulated rate; and R(t i , d) is the two-way instantaneous slant range. Forty bistatic SAR echo signals can be obtained for the BingSat-TOPIS system, and the radar data cube for a point target in ice can be expressed as (38)

A. Satellite Orbit
In Table I, the system parameters of BingSat-TOPIS are given. The Kepler Orbit can be explained by six Keplerian elements. There are three elements for the orientation of the orbital plane, which are inclination (i), argument of perigee (ω), and right ascension of the ascending node (Ω). The shape and size of ellipse can be defined by eccentricity (e), semimajor axis (a), and true anomaly (f ). These Keplerian elements determine the satellite position in a 3-D space.
Inclination i, semimajor axis a, and eccentricity e are fixed parameters, which argument of perigee ω, right ascension of the ascending node Ω, and true anomaly f varied with the time (t). The relationship between eccentric anomaly (E) and mean anomaly (M ) can be described as Mean anomaly M is relate to the epoch time, which is defined as where μ = 398601.2km 3 /s 2 is the Earth's gravitational constant, and a is the semimajor axis. Suppose that the satellite passes through the perigee at time t 0 = 0 so that Substituting (25) into (23), the eccentric anomaly E is obtained.
The relationship between the eccentric anomaly E and the true anomaly θ is given by In the orbital coordinate system, the satellite position vector can be described as where r is the distance between the satellite and the Earth's center, and r can be written as In orbital coordinate system, the satellite velocity vector, the first derivative of the position vector (with respect to t), can be  The first derivative of r , with respect to t, is given bẏ Based on the motion equation of satellite, the energy of particle per unit mass, h, can be written as Substituting (47) Using the coordinate transform equation, the vectors of the satellite position, velocity, and acceleration are described in the ECR coordinate system. The vectors of the satellite position, velocity, and acceleration in the ECR coordinate system are As we know, a circle is a special kind of ellipse (e = 0). The master satellite of BingSat-TOPIS is a circular orbit, and the slave satellites of BingSat-TOPIS are elliptical orbits.
The system parameters of BingSat-TOPIS are listed in Table I. By using these parameters, we could simulate the orbits of the master and 40 slave satellites.

B. Scene Setting
t c is the azimuth center time. According to Section IV-A, the vectors of the satellite position, velocity and acceleration can be calculated. Then, the master satellite's geocentric latitude and geocentric longitude of the subsatellite point at t c are given by where χ sm and φ sm are the master satellite's geocentric latitude and geocentric longitude, respectively.[x sm , y sm , z sm ] T is master satellite's position. The master satellite's geocentric longitude is the same as the master satellite's geodetic longitude, but its geodetic latitude is not the same as its geodetic latitude. The master satellite's geodetic latitude can be obtained where e c is the Earth's eccentricity, h is the ellipsoid height, and the radius of curvature in the prime vertical N T is given by Supposing that the scene's center is located under the Earth's surface 2000 m (h = −2000), whose geodetic latitude and geodetic longitude are the same as the subsatellite point at t c . Therefore, the position vector of scene's center in the ECR coordinate system is given by Many point targets could be placed in one scene, and each of them has the only position vector in the ECR coordinate system.
In order to make all the satellites just over the North Pole area, the azimuth time center t c = 79800s. Nine point targets are placed below the master satellite, whose geodetic latitude, geodetic longitude, and ellipsoid height are listed in Table II.
The spatial positions of nine point targets are shown in Fig. 7. The nine point targets are all ideal points, whose backscattering coefficients equal one.

C. Tomographic Imaging Steps and Results
In Fig. 8, the tomographic imaging processing steps are illustrated. As we know, there are 40 bistatic SAR echoes for BingSat-TOPIS. First, the 40 echoes should be range compressed by using the matched filter in frequency domain. For all the 40 range compressed data, the range interpolation operation should be performed. For an imaging area, it should be divided into many 2-D grids. Back projection (BP) algorithm [16] is used in track and we obtain the 2-D image, which includes calculation of the equivalent range history of bistatic SAR by a new model, phase compensation, and accumulation to every pixel. Then, 3-D grids are generated. Finally, BP is used in cross-track by using all the 40 2-D images and we obtain the 3-D imaging result. The BP algorithm can realize track and cross-track focusing without phase errors.
As shown in Fig. 9, the tomographic imaging results of nine point targets are displayed in two kinds of 3-D coordinate and one 2-D profile.
As we can see in Fig. 9(b), the coordinates of the nine points in track are all equal to zero. The other expression for a point is given by (cross-track, track/azimuth, and range). The point target located at scene center is (00, -2000). Compare to the scene center, the values of cross-track and track are given.
Three To examine the responses of nine point targets in more detail, the amplitudes of nine points' slices in track/azimuth, crosstrack, and range are shown in Figs. 10-12, respectively. All the three figures are shown after interpolation.
The resolutions of nine points are given in Table III. The peak sidelobe ratio (PSLR) and the integrated sidelobe ratio (ISLR) of nine points are given in Tables IV and V, respectively. All the results can meet our design indexes.

V. CONCLUSION
Tomographic imaging for orbital radar sounding of Earth's ice sheets is proposed in this article. Dozens of MirrorSAR Cube-Sats operating with SIMO mode could realize high range, track, and cross-track resolution. In this article, the ionospheric effects are ignored. In most cases (TEC value is less than 30TECU), image defocusing in track and cross-track can be ignored. Only during extreme ionospheric activity, it is necessary to consider ionospheric compensation by one value in track. By using the matched filter in the frequency domain, the range compressed data can be easily obtained. Then, the BP algorithm is used in both track and cross-track to get the tomographic image. We deduce the equivalent slant range of a new model, which can improve the computational efficiency by using it in the BP algorithm. Since 2018, she has been a Senior Engineer with the Qian Xuesen Laboratory of Space Technology, China Academy of Space Technology, Beijing. Her research interests include synthetic aperture radar (SAR) image processing, SAR image quality improvement, and SAR image geolocation without ground control points. Her current research interest include tomographic imaging for radar sounding. From 2014 to 2017, he held a postdoctoral position with Beihang University. He has been a Senior Engineer with the Qian Xuesen Laboratory of Space Technology, China Academy of Space Technology, Beijing. His research interests include innovative remote sensing systems and applications based on ultrawide swath imaging with synthetic aperture radar systems, orbital radar sounder, and the development of advanced signal algorithms.