Doppler Analysis of Forward Scattering Radar With Opportunistic Signals From LEO Satellites

Recently microwave communication signals from low Earth orbit (LEO) satellites have been proposed to be used to detect flying objects in a ground-based forward scattering passive radar setup. In this work the Doppler shift and the Doppler spread due to the flying object in such a system are derived directly from the Fresnel-Kirchhoff diffraction formula using a rotational Cartesian coordinate system to accommodate the continuously moving signal source from the LEO satellites. Coordinate transformations for flying objects including those in the atmosphere with rectilinear motion relative to the ground and those in space with circular orbits are investigated. Both analytical and simulation results demonstrate that the Doppler shift and Doppler spread are related to the velocity and volumetric profile of the flying object, thereby in future they may be used to help identify its size and velocity.


I. INTRODUCTION
Bi-static forward scattering radar has been widely investigated since it was first envisioned about seven decades ago [1]. By exploiting various opportunistic signals, such as GNSS signals [2], cosmic radio emissions [3], terrestrial broadcasting signals [4], DVB-S signals [5], recently ground based forward scattering (FS) passive radars have attracted a lot of attention for their capability in air object detection. Using opportunistic signals from satellites, CubeSat-based FS passive radars have also been proposed for space debris detection [6], [7].
Along with the deployment of Low Earth Orbit (LEO) satellite constellations for global Internet services such as Starlink, more and more Ku, Ka, and V band microwave signals will be available from the sky for communications. We have recently proposed to use such signals opportunistically to retrieve things in the sky, such as rainfall [8] and flying objects [9], with a ground-based FS passive radar setup. As pointed out in [8] and [9], advantages of such systems include low cost, availability of abundant opportunistic sources, and the networking capability.
The associate editor coordinating the review of this manuscript and approving it for publication was Vittorio Degli-Esposti . Doppler analysis is critical for detection and imaging with radars. The Doppler effects for FS radars have been analyzed in [10] by utilizing a two-path signal model with the assumption that the direct path signal can be removed by using a self-mixing heterodyne or an envelope detector [4], [11]. However, in a FS passive radar setup where opportunistic signals of a large bandwidth such as the communication signals are employed, a rigorous analysis has not been available. Furthermore, to date only static coordinate systems were used for research in FS bi-static radars [2], [3], [4], [5], [6], [10], [12]. When signals from LEO satellites are employed opportunistically, a static coordinate system cannot accommodate the continuously moving signal source, as in such a scenario the optical axis of the diffraction system is continuously moving.
In this work, we are focused on the Doppler effect of a flying object for ground based FS passive radars with opportunistic signals from LEO satellites, albeit the proposed method can also be applied to space-based radars and other signal sources. In contrast to previous work that used a simplified two-path signal propagation model [10], the Doppler shift due to a flying object is derived directly from the Fresnel-Kirchhoff diffraction model and Babinet's principle, thereby the assumption of the direct path signal removal is not VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ required, which is important for signals with large bandwidth thereby hard to separate the direct path from the indirect path, such as the communication signals of Starlink. Furthermore, from the Fresnel-Kirchhoff diffraction model, it will be shown that the Doppler spread (i.e. the spread of Doppler shifts. This happened when the size of the flying object is too large to be treated as one point) can be characterized by two bounds, which are related to the size and the velocity of the flying object. In contrast to the static coordinate system employed in the literature, in this work, to accommodate the continuously moving LEO satellite, we adopt a rotational coordinate system, which is attached to the link between the LEO satellite and the ground receiver, to make the optical axis of the Fresnel-Kirchhoff diffraction aligned with the link between the satellite transmitter and the ground receiver. As a result, the coordinate system that characterizes the motion of the flying object needs to be transformed into the rotational coordinate system. In this work, as an example, we will demonstrate coordinate transformations for two types of flying object motion: flying object in the atmosphere with rectilinear motion relative to the ground and flying object in space with circular orbit 1 .
To verify the proposed method, Ka band microwave communication signals from a LEO satellite propagate through a flying object are simulated. As the symbol rate of modern satellite communication systems are often much larger than the Doppler frequency, the amplitude of the received signal at ground station is averaged over a time window to eliminate the impact of the time variations of the data symbols. The analytical results of Doppler frequency and Doppler spread (i.e. the two bounds) are then compared with the spectrograms of the average amplitude of the received signals at the ground station. The impact of different window length, signal to noise power ratio (SNR) and satellite elevation angle is investigated through simulations. Simulation results show that the spectrograms of received signals are consistent with the analytical results for both flying object in the atmosphere in a rectilinear motion relative to the ground and that in space in a circular orbit. Particularly, the size and the velocity of the flying objects have an impact on the Doppler spectrum, thereby paving the way for the use of the FS passive radars with opportunistic signals from LEO satellites to detect relevant parameters of, for example, an airplane or a space debris, in future. This paper is organized as follows. The system model is introduced in Section II. The Doppler analysis is presented in Section III, followed by two examples of coordinate transformations in Section IV, one for fly object in the atmosphere, the other for that in space. Simulation results are presented in Section V, and conclusions are drawn in Section VI.

II. SYSTEM MODEL
As shown in Fig. 1, a rotational Cartesian coordinate system (x, y, z) is defined as follows. The ground receiver, which is fixed on the ground, is located at the origin. The x axis is parallel to the ground plane, and the y axis is along the link between the satellite and the ground receiver, which rotates around the ground receiver. Both the location of the ground receiver and the orbit trajectory of the LEO satellite is assumed to be known, 2 thereby the coordinate system (x, y, z) is assumed to be given. Assume that the distance between the satellite and the ground receiver is y S (t). Then The coordinate of the satellite is (0, y S (t), 0). (x , y , z ) is another rotational Cartesian coordinate system with the origin attached to the center of the flying object and the (x , y , z ) axes are parallel with the (x, y, z) axes, respectively.
Using the Fresnel-Kirchhoff diffraction formula and Babinet's principle, the complex amplitude of the received signal at the ground station is given by (1), as shown at the bottom of the next page, [14], [15]. In (1), A is the complex amplitude of the signal at the satellite transmitter, the wave number k = 2π/λ, where λ is the wavelength, the diffraction angles α 1 and α 2 are both approximated as 0, and where is the 2D image of the flying object projected onto the x − z plane. Note that here we assume that the thickness of the flying object in the y direction is negligible. Also note that (1) is a baseband equivalent model, thereby the Doppler shift of the carrier of the communication signals due to the movement of the LEO satellite only, which can be known a priori [16], is assumed to having been compensated and thereby be ignored in this work. Let the coordinate of the centroid of the flying object at time t under the (x, y, z) system be (x F (t), y F (t), z F (t)). Then the distance between a point (x , 0, z ) on the flying object and the ground receiver, and that between the point and the satellite are given by and respectively. Let Assume that r 1 (t) is much larger than x and z of interest. With the approximation of (2) can be approximated as follows.
where 2 The accuracy of the positioning of a LEO satellite can be as high as a few centimeters [13]. By the same token, let . Assume that r 2 (t) is much larger than x and z of interest, then (3) can be approximated as follows. where

III. DOPPLER ANALYSIS A. DOPPLER SHIFT
From (1), the phase difference for a point (x , z ) on the flying object to that of the direct path is given by With the approximation of r 1 by (4) and r 2 by (7), the Doppler shift due to the centroid of the flying object is given by (11), as shown at the bottom of the next page. As an example, with the same setup as in Section V for the simulation of the flying object in the atmosphere, the Doppler shifts computed by (11) at a LEO satellite elevation angle of 65 degree are shown in Fig. 2 for the speeds of the flying object being 300m/s, 150m/s, 0m/s, −150m/s and −300m/s. It can be seen that the Doppler shifts change linearly in time and the slopes are related to the speed of the flying object. Particularly, a speed of 150m/s leads to an almost zero Doppler shift, as the relative velocity (observed at the flying object) between the flying object and the LEO satellite is small. Let . With the use of (5), (6), (8), (9), and (11), the Doppler shift due to a point (x , z ) on the flying object is given by (12), where dη(t) dt can be computed using (13).
In general for a flying object in the atmosphere, r 2 (t) is at least one order larger than r 1 (t), thereby r −3 2 (t) is three order smaller than r −3 1 (t). As a result, the second term of (13) can be ignored, as it is proportional to r −3 2 (t) while the first term is proportional to r −3 1 (t). 3 In similar token, the third term in (12) is much smaller than the second term. As a result, f x ,z (t) can then be approximated as follows. 4 3 From (13), , which, particularly for relative high satellite elevation angles as in broadband LEO satellite systems, is one or two order larger than However, the order of r −3 2 (t) and r −3 1 (t) dominates. 4 If the flying object is in space (i.e. in orbit similar to a LEO satellite or a space debris), then the approximation of (14) is still valid if the flying object and the LEO satellite fly in the same direction. If they fly in the opposite direction, then the apporoximation of (14) is in general not valid. In such a scenario, however, the time window for observing the flying object will be too short to make successful estimation of the Doppler shift or Doppler spread possible. As an example, with the same setup as the simulations for Fig. 2, Fig. 3 shows the Doppler shifts of microwave signals from a LEO satellite due to the point at one of the four corners of a rectangular flying object of size {a = 20m, b = 40m} (a and b are defined in Section V) in the atmosphere at an elevation angle of 65 degree, which demonstrate that the difference between (12) and the approximation of (14) is negligible.

B. BOUNDS TO CHARACTERIZE DOPPLER SPREAD
For a digital satellite communication system such as Starlink, the complex amplitude of the transmitted signal by the satellite is not a constant, thereby A in (1) needs to be replaced with a time-varying complex-valued variable to accommodate the different data symbols in time that are sent by the transmitter. 5 Let the first term of (1) be A(t)e jϕ(t) . Then the second term of (1) can be represented as A(t)e jϕ(t) A x ,z (t)e jϕ x ,z (t) dx dz , where A x ,z (t) and ϕ x ,z (t) are the amplitude and the phase associated with point (x , z ) on the flying object. As a result, the RF signal can be represented as (15), shown at the bottom 5 Note that this does not impact the Doppler shift analysis in Section III-A, as the phase difference given by (10) is still valid. of the next page, where f c is the carrier frequency. Note that in general the time scale for the change of A(t) and ϕ(t) (for a satellite communication system, it is in the order of microsecond or less) is much smaller than that for the change of A x ,z (t) and ϕ x ,z (t) (in the order of millisecond).
Let the amplitude of S(t) in (15) be S A (t). Then from (16), as shown at the bottom of the next page, and Appendix A, we have As A(t) changes much faster than A x ,z (t) and ϕ x ,z (t), the time-variation of A(t) can be averaged out over a time window sufficiently large using, e.g., a low pass filter. Without loss of generality, assume that E(A(t)) = A, where E(.) is timeaveraging of (.). Then From (18), it can be seen that after time-averaging of the amplitude of the received signal, positive frequency components and negative frequency components cannot be differentiated. Instead, |f d (t)| and |f x ,z (t)| can be measured using, e.g. the spectrogram of E(S A (t)). Based on (12) or (14), different point in a flying object may produce different Doppler shift, thereby a flying object of adequate size will produce a spread of Doppler shifts. We use the following two bounds to characterize the Doppler spread: and As an example, Fig. 4 shows the Doppler shift, bound1 and bound2 of a rectangular flying object of size {a = 20m, b = 40m} in the atmosphere at an elevation angle of 65 degree. The speeds of the flying object are 300m/s, −150m/s, and −300m/s. It can be seen that, given an elevation angle, size and speed of the flying object, the difference between bound1 and bound2 in general does not change with time. Therefore, the difference between bound1 and bound2 is defined as the Doppler spread as follows.  Fig. 5. It can be seen that the Doppler spread is related to both the size and the speed of the flying object. The larger the size of the flying object, the larger the Doppler spread. The Doppler spread is minimum if the flying object is at a speed of about 150m/s (at this speed for a satellite elevation angle of 65 degree, the relative movement between the satellite and the fling object is small), and increases along with the increase of the speed difference with 150m/s.

IV. COORDINATE TRANSFORMATIONS
The variables in (11) and (14) are defined under the rotational coordinate (x, y, z) system. However, in practice, the coordinates and the velocities of a flying object are often defined under a different coordinate, e.g. the trajectory of an airplane is often defined under a coordinate that is fixed on the ground, and the trajectory of a satellite or a space debris is often defined under the Earth-centered, Earth-fixed coordinate system. In the following, how to convert the velocities under a coordinate system different from the rotational coordinate (x, y, z) system into v x F , v y F , v z F (in the following, v x , v y , v z will be used for notational simplicity) and v y S in (11) and (14) will be demonstrated using two examples.

A. FLYING OBJECT IN THE ATMOSPHERE
As shown in Fig. 6, let (x, y G , z G ) be a Cartesian coordinate system fixed on the ground with the x − y G plane being the  . Relationship between coordinate system fixed on the ground (x, y G , z G ) and rotational coordinate system (x, y , z). Both systems have the same x axis, which is perpendicular to the yOz (or y G Oz G ) plane.
ground plane and the x axis and the origin being the same as those in the (x, y, z) coordinate system. Let θ be the angle between axis y in the (x, y, z) system and axis y G in the (x, y G , z G ) system. Then a point (y G , z G ) in (x, y G , z G ) can be converted into a point (y, z) in the (x, y, z) system as follows (note that the x coordinates are the same): Let v y = dy dt , v z = dz dt , v y G = dy G dt and v z G = dz G dt . Then from (21), v y and v z can be computed from y G , z G , v y G , v z G , θ and To find dθ dt , let us consider a circular LEO orbit with Earth's rotation being ignored [16]. The orbit height is assumed to be h Sat . As shown in Fig. 7, the relationship between the elevation angle θ of a satellite relative to a ground station and the angle α relative to Earth center is as follows.
where r E = 6.371 × 10 6 m is the radius of the Earth. As a result, dθ dt can be computed by (24), as shown at the bottom of the next page. In (24), dα dt = 2π T Sat and the period of the satellite T Sat = 1 2π (h Sat +r E ) 3 GM Earth , where G = 6.6743e−11m 3 kg −1 s −2 is the universal gravitational constant and M Earth = 5.9722e24kg is the mass of the Earth.
To use (11), v y S (t) also needs to be found. From Fig. 7, it can be found that v y S (t) can be computed by (25), as shown at the bottom of the next page.

B. FLYING OBJECT IN SPACE
For simplicity, we assume that the flying object of interest (e.g. a space debris) and the LEO satellite are orbiting on the same plane but with different orbit height. As shown in Fig. 8, the coordinates of the debris under the rotational coordinate system (x, y, z) are given as follows: where d debris is the distance between the ground receiver and the debris, θ and θ debris are the elevation angles of the LEO satellite and the debris, respectively. As a result, v x (t) = 0. v y (t) and v z (t) are given by (26) and (27), as shown at the bottom of the page, respectively. In (26) and (27), dθ debris (t) dt and dθ(t) dt can be computed using (24); d debris (t) dt can be computed using (25).
In Section V, we assume that at time 0, the satellite and the space debris have the same elevation angle θ . Given θ at time 0, α can be obtained as follows. If θ ≤ π/2, then α is given by (28), as shown at the bottom of the page. Otherwise, if θ > π/2, then α is given by (29), as shown at the bottom of the page.
Given the α at time 0, α at time t can then be obtained as α(t) = α(t = 0) + 2π t/T , where T is the period of the LEO satellite or the debris. θ at time t can then be obtained from α(t) using (23).

V. SIMULATION RESULTS
In practice Doppler shifts are often measured using the spectrograms, e.g. as in [4]. In this work, to verify the analytical results of (11), (19) and (20), which are only related to the coordinates and velocities of the flying object and the LEO satellite, we use the spectrograms of the amplitude of the simulated received microwave signals from a LEO satellite.
In the simulations, we assume a LEO satellite with a circular orbit of height h Sat = 550km. When passing over the ground receiver, the LEO satellite is assumed to be vertically above it. The wavelength of the microwave signals is assumed to be λ = 1 cm (i.e. f c = 30GHz). For simplicity the flying (29) VOLUME 10, 2022 object is modeled by a rectangular slab with width a meters and length b meters.
The simulations are carried out as follows. A constant signal or a communication signal is generated from the LEO satellite. The complex amplitude of the received signal at the ground receiver is computed with (1) considering the shape and the velocity of the flying object and the orbit of the LEO satellite. The amplitude of the received signal is computed from the complex amplitude of the received signals (additive white noise is added when necessary), which is then analyzed with the short Fast Fourier Transform to produce the spectrogram. The results of the spectrogram is then compared with the analytical results of (11), (19) and (20). We first simulate a flying object in the atmosphere. As shown in Fig. 6, it flies on the y G Oz G plane with a constant height h = 10,000m above the ground and a translational movement with a constant speed along the y G direction. The flying object is assumed be always parallel with the ground, and its width is in the x direction and its length in the y G direction. As a result, at time t, the length of the image of the flying object projected onto the xOz plane in the z direction is bsin(θ (t)). The flying object is assumed to be across the link between the satellite and the ground receiver at t = 0. A flying object of {a = 10m, b = 20m} with a speed of 250m/s (the flying object and the satellite move in the same direction) and −250m/s (the flying object and the satellite move in the opposite direction) are simulated to obtain the amplitude of the received signal. The complex amplitude of the transmitted signal at the LEO satellite is assumed to be a constant in time. It can be seen from Fig. 9 that the Doppler shift and Doppler spread derived in this work are consistent with the spectrograms of the amplitude of the received signals at the ground receiver (i.e., spectrograms of |U (t)|), which are related to speed of the flying object and the satellite elevation angles. When the flying object is across the communication link, the Doppler frequency is about zero Hertz, and then increases with time linearly. Specifically, the lower the satellite elevation angle, the smaller the slope of the increase. In the time domain, the average amplitude of the received signal increases long with the increase of the satellite elevation angle, as the distance between the LEO satellite and the ground receiver becomes shorter, leading to the decrease of the free space path loss.
We then consider communication signals with additive white Gaussian noise (AWGN). A communication signal with standard 16QAM modulation scheme [17] is assumed. To accommodate the noise, the received signal is now given by R(t) = U (t) + N (t), where N (t) is a complex valued random process to model the white Gaussian noise. At the ground receiver, R(t) is assumed to be sampled every data symbol duration. At different sample time t n , the noise N (t n ) are independent and identically distributed Gaussian random variables. The amplitude of received signals is then approximated as A ≈ W n=1 |R(t n )|, where W is the window length. For an SNR 6 of 0 dB, 10 dB and 20 dB and a window length W = 100 and W = 1000, Fig. 10 shows the spectrograms for a flying object of {a = 20m, b = 40m} with a speed of 250m/s for a satellite elevation angle of θ(t = 0) = 65 • . While spectrograms can be used to characterize the Doppler effects, it can be seen from Fig. 10 that both the window size and the SNR impact the performance of the spectrograms, which can be improved by increasing the window size (e.g. comparing Fig. 10a and Fig. 10c) and/or increasing SNR (e.g. comparing Fig. 10c and Fig. 10f). Particularly, for low SNR and small window size (e.g. Fig. 10a), the Doppler effects are barely visible through spectrograms. Overall, it can be seen that the spectrograms are consistent with the Doppler shift and the Doppler spread (i.e. bound1 and bound2) derived in this work. It can also be seen that the Doppler effect becomes more prominent by increasing either window length or SNR. Compared Fig. 9e with Fig. 10f, it can be seen that increasing the size of the flying object will make Doppler spread larger, which is consistent with the results of Fig. 5.
Finally we simulate a rectangular flying object (space debris) in orbit. It is assumed that the LEO satellite and the space debris are on the same plane, flying in the same direction, but with different orbit height. Fig. 11 shows the spectrograms for a space debris with a = 20m (in the x direction), b = 40m (on the yOz plane in the direction that is perpendicular to the link between the space debris and the ground receiver), height of 500km and 540km, with SNR = 10 dB, W = 1000, and θ (t = 0) = 40 o , 65 o , 90 o . It can be seen that the spectrograms are consistent with the analytical results of Doppler shift and Doppler spread. It can also be seen that the closer of the orbit of the object to the LEO satellite, the more prominent the spectrograms, but the Doppler shift and Doppler spread becomes smaller, as the relative speed between the LEO satellite and the space debris becomes smaller. Same as flying object in the atmosphere, for relative low satellite elevation angle the Doppler frequencies are smaller and in the time domain the average amplitude of the received signal is smaller, due to the larger free space path loss.

VI. CONCLUSION
The Doppler shift and the Doppler spread due to a flying object for microwave signals from LEO satellites in a bistatic forward scattering passive radar setup has been derived in this work with the use of a rotational coordinate that is attached to the communication link between the LEO satellite and the ground receiver. Consistent with the Doppler analysis, simulation results of the spectrograms of the amplitude of the received signals for a flying object in the atmosphere and a space debris in space have verified that the Doppler spread and the Doppler shift are related to the size and velocity of the flying object, thereby they may be used to detect relevant parameters of the flying object in future.

APPENDIX A
From (15), S A (t) can be derived as (16). Note that in (16), for the third line, the following inequality is applied.
For the fourth line, the approximation of √ 1 + x ≈ 1 + x/2 for small x is applied. For the fifth line, the second order term 1 2 A 2 x ,z (t)dx dz is ignored, as in general it is much smaller than the absolute value of the first order term.
From the second line of (16), it can be seen that S A (t) ≥ A(t) 1 + A x ,z (t)cos(ϕ x ,z (t))dx dz ) .