A Novel Ship Imaging Method with Multiple Sinusoidal Functions to Match Rotation E ﬀ ects in Geosynchronous SAR

: Geosynchronous Synthetic Aperture Radar (GEO SAR) has a very long Coherent Processing Interval (in the order of hundreds of seconds) compared with other SAR platforms. Thus, the current methods of rotation e ﬀ ect matching and ship imaging that operate within a relatively short Coherent Processing Interval (in the order of seconds) are obviously not applicable. To address this problem, a novel ship imaging method with multiple sinusoidal functions matching for rotation e ﬀ ects is proposed for GEO SAR. Firstly, the inﬂuence of the rotational motion of a ship on the slant range is analyzed. It can be matched with the sum of multiple sinusoidal functions, and the signal model of a ship with rotational motion is given. Then, multiple sinusoidal functions for the matching-based ship imaging method are proposed, and their procedures are presented as follows: (1) The Generalized Keystone Transform and Generalized Dechirp Process (GKTGDP) is modiﬁed to compensate for the range migration and phase caused by the motion of GEO SAR. Then, the signal is focused at the frequencies of sinusoidal functions, and the frequencies can be matched. (2) From the matched frequencies, the other parameters of sinusoidal functions can be matched by parameter searching. (3) Based on the matched results, the Back Projection Algorithm (BPA) is used to take an image of the ship with rotational motion. Finally, the e ﬀ ectiveness of the proposed method is veriﬁed by numerical experiments.


Introduction
Synthetic Aperture Radar (SAR) can work all day and in all weather conditions to obtain high-resolution remote sensing images [1]. It has become the main approach to ship surveillance [2]. Studies on ship surveillance have mainly focused on ship target detection and estimation [3], multichannel sea clutter suppression [4], and target refocusing [5]. However, currently, almost all SAR satellites run on the Low Earth Orbit (LEO) with an orbit height of around 500-1000 km. Limited by their orbit, an LEO SAR has a coarse time resolution (days to weeks) and limited swath coverage (tens of kilometers).
Geosynchronous SAR (GEO SAR) is an SAR satellite with an orbital height of around 36,500 km that can provide daily coverage for approximately one-third of the Earth's surface [6,7]. Therefore, Meanwhile, there is a curved synthetic aperture trajectory and a non-negligible signal round-trip delay in GEO SAR. Thus, the existing GKTGDP cannot match the rotational motion of a ship in GEO SAR.
To realize rotation effect matching and ship imaging in GEO SAR, a novel ship imaging approach using multiple sinusoidal functions to match the rotation effect is proposed in this paper. In Section 2, the influence of rotation on the slant range is analyzed and can be matched by the sum of three sinusoidal functions. Each of the sinusoidal functions has an unknown frequency, amplitude, and initial phase. Then, the signal model of a rotating ship in GEO SAR is given. In Section 3, the multiple sinusoidal function matching-based rotational ship method is proposed, whose procedures can be presented as follows: (1) The frequencies of sinusoidal functions are matched. The GKT is modified to compensate for the range migration caused by the motion of the GEO platform. Meanwhile, the curved trajectory and non-negligible signal round-trip delay are compensated for. Then, the GDP is used to focus the signal at the frequencies of sinusoidal functions, and the frequencies can be matched, i.e., the GKTGDP is modified to realize frequency matching. (2) The amplitude and initial phase of sinusoidal functions are matched. The phase of a signal is compensated for by searching for the amplitude and initial phase. When the signal is completely focused at the Doppler center frequency, amplitude and initial phase matching are realized. (3) The imaging method of a rotating ship is carried out. Based on the matched results of sinusoidal functions, the Back Projection Algorithm (BPA) is used to take an image of the ship with rotational motion. Most of the existing rotational ship imaging algorithms are based on the assumption that the translational motion has been compensated for and each degree of rotation is approximately a sinusoidal function [35], and this paper has the same assumption. Comparatively, the CPI of GEO SAR with low orbit inclination is much longer than that of GEO SAR with high orbit inclination (in the order of dozens of minutes) [9,10], and the rotation is hardly approximated as a sinusoidal function. Therefore, the proposed method is, currently, difficult to apply in GEO SAR with low orbit inclination. In Section 4, simulations are carried out to verify the proposed method. Section 5 concludes the paper.

The Geometry of a Rotating Ship
The rotational effect on the slant range of a scatter point of a ship is analyzed in this section. Firstly, the rotational motion and position of a scatter point are proposed for the ship's coordinate system. Then, the influence of rotational motion on a scatter point is transformed to the scene's coordinate system. Finally, the influence of rotational motion is projected to the slant range of the SAR. Figure 1 shows the general geometry of a rotating ship in GEO SAR. The ship's fulcrum is usually made up of the gravity center O and the origin of the ship coordinate system OBPZ. The geometry of the scatter point of a ship is expressed by this coordinate system. In OBPZ, the positive Band P-axes point to the forward direction and the port side of a ship, respectively. The Z-axis is the ship vertical axis. It is assumed that the ship only has 3D rotational motion, defined by roll, pitch, and yaw, which is measured counterclockwise from the B-, P-, and Z-axes, respectively. Based on the ship dynamics, the main forcing function driving ship rotational motions is waves at the sea surface. These waves interact with a ship's buoyancy and hydrodynamic forces to induce oscillatory motions whose angle tends to be sinusoidal [47,48] and can be expressed as a sinusoidal function [36]: θ r,p,y (t) = Θ r,p,y sin(2π f r,p,y t + ϕ r,p,y ) (1) where θ r,p,y represents the time-varying angles of 3D rotation such as roll, pitch, and yaw; Θ r,p,y , f r,p,y , and ϕ r,p,y are the amplitudes, frequencies, and initial phase of 3D rotation, respectively; and t is the slow time. Each scatter point has a coordinate at the aperture center moment (ACM). Due to the rotational motion of the ship, the coordinates change with t. To analyze the change in a scatter point's coordinate with t, based on Euler's rotation theorem, the transformation matrices of the roll, pitch, and yaw angles are shown as where θ r , θ p , and θ y represent the roll angle θ r (t), pitch angle θ p (t), and yaw angle θ y (t) at t, respectively, whose changes follow (1).
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 27 point's coordinate with t , based on Euler's rotation theorem, the transformation matrices of the roll, pitch, and yaw angles are shown as where θ r , θ p , and θ y represent the roll angle ( ) θ r t , pitch angle ( ) θ p t , and yaw angle ( ) θ y t at t , respectively, whose changes follow (1). OXYZ is the scene coordinate system, in which the Y -axis represents the direction of the SAR platform's velocity at the ACM. The X -axis is perpendicular to YOZ and is determined by the right-hand rule. Then, the direction angle ϕ c is the angle between the B -and X -axes. The transformation matrix from OBPZ to OXYZ can be expressed as OX Y Z is the imaging coordinate system, where the ′ X -axis is the direction of the slant range of SAR; the ′ Y -axis is perpendicular to the slant range and belongs to XOY ; the ′ Z -axis is perpendicular to ′ ′ X OY and is decided by the right-hand rule. The position of the scatter point projected to the slant range and the azimuth direction at the ACM in the coordinate system ′ ′ ′ OX Y Z OXYZ is the scene coordinate system, in which the Y-axis represents the direction of the SAR platform's velocity at the ACM. The X-axis is perpendicular to YOZ and is determined by the right-hand rule. Then, the direction angle ϕ c is the angle between the Band X-axes. The transformation matrix from OBPZ to OXYZ can be expressed as OX Y Z is the imaging coordinate system, where the X -axis is the direction of the slant range of SAR; the Y -axis is perpendicular to the slant range and belongs to XOY; the Z -axis is perpendicular to X OY and is decided by the right-hand rule. The position of the scatter point projected to the slant range and the azimuth direction at the ACM in the coordinate system OX Y Z . The transformation Remote Sens. 2020, 12, 2249 5 of 24 matrix from OXYZ to OX Y Z is a product of the transformation matrix of the grazing angle and the angle between the X-axis and the XOY plane, and can be written as cos ϕ g cos ϕ s cos ϕ g sin ϕ s − sin ϕ g − sin ϕ s cos ϕ s 0 sin ϕ g cos ϕ s sin ϕ g sin ϕ s cos ϕ g where angle ϕ g is the grazing angle; angle ϕ s denotes the angle between the X-axis and the projection of the X -axis on the XOY plane in the ACM. From (2) to (6), the coordinates for the ith scatter point in the OX Y Z coordinate system can be expressed as where b i p i z i T are the ith scatter point coordinates of a ship in the OBPZ; y i (t) and z i (t) are the projections of the ith scatter point on the Y -and Z -axes at t, respectively; and r i (t) is the projection of the ith scatter point on the slant range (X -axis) at t and can be rewritten as follows: r i (t) = cos ϕ g cos ϕ cs b i cos θ p cos θ y + p i cos θ p sin θ y + z i sin θ p where angle ϕ cs = ϕ c + ϕ s is the angle between the slant range and the B-axis. From (8), it can be seen that r i (t) is influenced by the 3D rotational motion, and it is difficult to match the rotation effect directly.

The Influence of the Rotation Effect
From (2) to (4), each degree of freedom of rotational motion can be expressed as a nested sinusoidal function: sin Θ r,p,y sin(2π f r,p,y t + ϕ r,p,y ) cos Θ r,p,y sin(2π f r,p,y t + ϕ r,p,y ) (9) where mathematical symbols have the same meaning as described in (1). Based on Jacobi-Anger expansion, which is the Bessel function of the first kind, (9) can be rewritten as sin Θ r,p,y sin(2π f r,p,y t + ϕ r,p,y ) cos Θ r,p,y sin(2π f r,p,y t + ϕ r,p,y ) where J α (Θ r,p,y ) is the Bessel function of α order and can be expressed as As shown in (8), the influence of rotational motion can be expressed as the sum of the products of the sinusoidal and co-sinusoidal functions shown in (10). Fortunately, the products of the sinusoidal Remote Sens. 2020, 12, 2249 6 of 24 and co-sinusoidal functions can be rewritten as the sum of them. Therefore, (8) can be rewritten as the sum of sinusoidal functions: where f n is the frequency of the nth sinusoidal function; ϕ n is the initial phase of the nth sinusoidal function; n is the number of the sinusoidal function; and N p is the total number of the sinusoidal functions. From (8), it can be seen that the positions of different scatter points are not the same, i.e., b i , p i , and z i are not the same between scatter points. Therefore, the amplitudes of the sinusoidal functions of different scatter points are different, and Θ i,n represents the amplitude of the nth sinusoidal function of the ith scatter point. The influence of rotation can be matched by (12), and the residual error of the matching result is minimal when N p = 3 in (12) (see Appendix A for the details). It can be seen that when the radar illuminates the ship, three rotational motions that trend to sinusoidal have been projected into the slant range. Although each rotation has a change in its scale after projection, there are still three sinusoidal functions projected into the slant range. Therefore, the rotation effect of a ship can be matched with the sum of three sinusoidal functions.
It is well known that the smaller the phase error is, the better the focusing effect is. Therefore, the rotational motion of a ship can be matched by gradually increasing the number of sinusoidal functions. When the total number is equal to three, the residual phase error after matching is minimal. Figure 2 shows the geometry of a ship with rotational motion in GEO SAR. GEO SAR operates in a Geosynchronous orbit and works in the side-looking mode. The X -axis is the direction from the GEO platform to the scatter point of a ship; the Y -axis is the azimuth direction; and the Z -axis is perpendicular to X OY plane. From (12), the slant range of a scatter point at t can be expressed as

The Signal Model of a Ship with Rotational Motion
where r O (t) is the slant range from GEO SAR to a scatter point of a ship, which is expressed as a fourth-order Taylor series and compensates for the signal round-trip delay [25]. It can be obtained as where r 0 , k 1 , k 2 , k 3 , and k 4 are the zero-order to fourth-order coefficients of the Taylor series, respectively; ∆r 0 , ∆k 1 , ∆k 2 and ∆k 3 are the zero-order to third-order coefficients of the Taylor series for the relative motion within the signal round-trip delay. In [25], the motion parameters of a ship with translational motion can be estimated in GEO SAR. Therefore, we assume that the translational motion has been compensated for and the coefficients of r O (t) are known from the motion of the GEO platform in this paper. As shown in (12) and Section 2.2, the rotation effect of a ship can be matched by the sum of three sinusoidal functions, and (13) can be rewritten as The influence of rotational motion on the range cell migration is analyzed in Appendix B. From the analysis result and the optimal resolution of GEO SAR (ρ r_best = 5m), it is shown that the rotational motion is less than half of the optimal resolution, and its influence does not cause range cell migration in GEO SAR.
Remote Sens. 2020, 12, 2249 7 of 24 respectively; Δ 0 r , Δ 1 k , Δ 2 k and Δ 3 k are the zero-order to third-order coefficients of the Taylor series for the relative motion within the signal round-trip delay. In [25], the motion parameters of a ship with translational motion can be estimated in GEO SAR. Therefore, we assume that the translational motion has been compensated for and the coefficients of ( ) O r t are known from the motion of the GEO platform in this paper. As shown in (12) and Section 2.2, the rotation effect of a ship can be matched by the sum of three sinusoidal functions, and (13) can be rewritten as The influence of rotational motion on the range cell migration is analyzed in Appendix B. From the analysis result and the optimal resolution of GEO SAR ( ρ = _ 5 r best m ), it is shown that the rotational motion is less than half of the optimal resolution, and its influence does not cause range cell migration in GEO SAR.
So, after the range compression, the echo signal model can be given as where γ  So, after the range compression, the echo signal model can be given as where B s = γT p , T p , and γ are the pulse bandwidth, the pulse duration of the transmitting signal, and the frequency modulation rate of the linear frequency modulated (LFM) signal, respectively; A rm = T p B s A r is the amplitude after the range compression and A r is the amplitude before the range compression.

The Rotation Effect Matching and Ship Imaging Procedure
As shown in Section 2, the influence of rotational motion can be matched by the sum of three sinusoidal functions. The multiple sinusoidal function matching-based rotational ship imaging method is summarized in Figure 3. It comprises the following three steps: Step 1: Match the frequencies of sinusoidal function After range compression, the echo signal is as shown in (16). Firstly, the range migration caused by the motion of the GEO platform is corrected with the modified GKT (see Section 3.2). Then, the phase caused by the platform motion is compensated for with the GDP. Because the rotational motion can be matched by the sum of three sinusoidal functions, the signal is focused at the frequencies of three sinusoidal functions in the Range-Doppler (RD) domain, and the frequenciesf n , n = 1, 2, 3 can be matched by the position where the signal has focused.
Step 2: Match the amplitudes and initial phases of sinusoidal function Although the frequencies of sinusoidal functionf n are matched in Step 1, the amplitudes Θ i,n and initial phases ϕ n still need to be matched. Because the intensities of signal atf n are different, the sinusoidal functions with the maximal and minimal intensity are named as the 1st and 3rd sinusoidal functions, respectively. The rest of sinusoidal function is named as the 2nd sinusoidal function. Firstly, Θ i,1 and ϕ 1 of the 1st sinusoidal function, are searched, and the searching values are used to compensate for the influence of the 1st sinusoidal function. For different searching values, the intensities of compensation results atf 1 are different. Since the intensity atf 1 is minimal, the influence of the 1st sinusoidal function is matched, and the searching values of the amplitude and initial phase are the matched values ofΘ i,1 andφ 1 , respectively. Then, the same approach is taken to compensate for the influences of the 2nd and 3rd sinusoidal functions. Since the influences of all three sinusoidal

Matching the Frequencies of Sinusoidal Function
The Keystone Transform (KT) is a common method to eliminate the effects of linear range migration [49]. Furthermore, the Generalized KT (GKT) has been proposed to eliminate high-order range migration [46]. In this Section, the GKT is modified to correct for the range migration caused by the motion of the GEO platform. The GDP is used to compensate for the phase caused by the platform's motion, and the frequencies of sinusoidal function can be matched in the RD domain. After transforming (16) into the slow time-range frequency domain, the signal can be expressed as where τ f is the range frequency corresponding to the fast time τ ; c f is the carrier frequency; and 1 A is the amplitude of the signal in the slow time-range frequency domain. Substituting (14) and (15) into (17) yields However, because different scatter points have the same frequencies and phases but different sinusoidal function amplitudes Θ i,n , the amplitudes of other scatter points need to be matched. Similarly, to match the sinusoidal functions of the ith scatter point, the Θ i,1 of the 1st sinusoidal function is searched in the range cell where the scatter point is, and the searching value is used to compensate for the influence of the 1st sinusoidal function. With the different searching value, the intensities of compensation results atf 1 are different. Since the intensity atf 1 is minimal, the influence of the 1st sinusoidal function is compensated for, and the searching value of the amplitude is the matched value ofΘ i,1 . Then, the Θ i,2 of the 2nd sinusoidal function is searched in the range cell. Since the intensity at f 2 is minimal, the searching value of the amplitude is the matched with the value ofΘ i,2 . Then, Θ i,3 was searched in the same way. Since the influences of all three sinusoidal functions are compensated for, the amplitudeΘ i,n of the scatter point is matched. Furthermore, when the amplitudes of all scatter points are matched, the influence of rotational motion matching is realized.
Note that scatter points' focusing in Steps 1 and 2 are in the RD domain and to match the parameters of sinusoidal functions. To realize the imaging of a ship, an imaging algorithm needs to focus each scatter point with the matched parameters in a 2D time domain.
Step 3: Ship imaging with the BPA In Steps 1 and 2, the f n , ϕ n , and Θ i,n of each scatter point of a ship are matched. Based on the matched values, a ship can be imaged by BPA.

Matching the Frequencies of Sinusoidal Function
The Keystone Transform (KT) is a common method to eliminate the effects of linear range migration [49]. Furthermore, the Generalized KT (GKT) has been proposed to eliminate high-order range migration [46]. In this Section, the GKT is modified to correct for the range migration caused by the motion of the GEO platform. The GDP is used to compensate for the phase caused by the Remote Sens. 2020, 12, 2249 9 of 24 platform's motion, and the frequencies of sinusoidal function can be matched in the RD domain. After transforming (16) into the slow time-range frequency domain, the signal can be expressed as where f τ is the range frequency corresponding to the fast time τ; f c is the carrier frequency; and A 1 is the amplitude of the signal in the slow time-range frequency domain. Substituting (14) and (15) into (17) yields where To correct the fourth-order range migration caused by k 4 , the fourth-order KT is performed along the slow time domain and the scaling formula can be written as Substituting (19) into (18) yields where λ = c/ f c . Correspondingly, the fourth-order phase can be compensated by following function By multiplying (20) with (21), we have Similarly, the third-order range migration and phase caused by (k 3 + ∆k 3 ) are compensated by the third-order KT and phase compensation, respectively. The third-order scaling formula for KT is given in (23); the third-order phase compensation function is given in (24): After the third-order KT and phase compensation, (22) can be rewritten as Then, to compensate for the range migration and phase caused by (k 2 + ∆k 2 ), the second-order scaling formula for KT and phase compensation function are expressed as (26) and (27), respectively: Because the motion parameters of GEO SAR are already known, the range migration and phase caused by (k 1 + ∆k 1 ) can be compensated for. The first-order scaling formula for KT and the phase compensation function are expressed as (28) and (29), respectively: After second-and first-order range migration and phase are compensated for, (25) can be rewritten as After performing an inverse fast Fourier transform (IFFT) along the range frequency domain, (30) can be expressed as where A 4 = A 3 sin c B s τ − 2 r 0 + ∆r 0 /c exp −j4π r 0 + ∆r 0 /λ ; A 3 is the amplitude of the signal after performing an IFFT. From (12), the influence of rotational motion can be matched by the sum of three sinusoidal functions, and (31) can be rewritten as Using (18) to (32), the first-to fourth-order range migration and phase are compensated for by KT and the phase compensation function, respectively. Each order of KT is collectively referred to as the generalized KT (GKT). Correspondingly, the first-to fourth-order phases need to be compensated for by the phase compensation function, which is similar to the Dechirp process (DP). Therefore, each order of phase compensation is collectively referred to as the generalized DP (GDP). Furthermore, each sinusoidal function can be expressed as a Bessel series [41].
where m is the number of harmonics, and α and β represent any amplitude and phase, respectively. From (33), (32) can be rewritten as Then, after transforming (34) into the RD domain, the signal will be focused at the frequency m f n [50][51][52]. A schematic diagram is shown in Figure 4. After GKT and GDP, the signal can be expressed as shown in (34) and is focused at the frequencies of sinusoidal function and the integral multiples of frequencies in the RD domain. Furthermore, from the Bessel function, the signal at ± f n is a fundamental component, and its intensity is greater. Meanwhile, the signal at m f n ,m = ±2, ±3, · · · are harmonic components, and their intensities are weak. Therefore, because the harmonic components are weak, the frequencies of sinusoidal function ± f n are easy to distinguish and can be matched by the position where the fundamental component of signal is focused, and the matched values of frequency aref n , n = 1, 2, 3.

Matching the Amplitudes and Initial Phases of Sinusoidal Function
From Section 3.2, after GKT and GDP, the signal can be focused on the frequencies of sinusoidal function, and the frequencies can be matched. However, from (15), the amplitudes Θ , i n and initial phases ϕ n of sinusoidal function also need to be matched.
Obviously, each compensation of the sinusoidal function can make the signal better focused at the Doppler center frequency dc f . A schematic diagram is shown in Figure 5.

Matching the Amplitudes and Initial Phases of Sinusoidal Function
From Section 3.2, after GKT and GDP, the signal can be focused on the frequencies of sinusoidal function, and the frequencies can be matched. However, from (15), the amplitudes Θ i,n and initial phases ϕ n of sinusoidal function also need to be matched.
Obviously, each compensation of the sinusoidal function can make the signal better focused at the Doppler center frequency f dc . A schematic diagram is shown in Figure 5.  Based on (34), searching parameters and phase compensation with searching values for each sinusoidal function in a range cell are shown in Figure 5. Firstly, the 1st sinusoidal function of (32) is compensated for by searching Θ i,1 and ϕ 1 . When the intensity atf 1 is minimal, the influence of the 1st sinusoidal function is compensated for and the searching values of the amplitude and initial phase are the matched values ofΘ i,1 andφ 1 , respectively. Then, the 2nd and 3rd sinusoidal functions are compensated for with the same method. The compensation function for the influence of the sinusoidal function can be expressed as M s,n t 1 = exp j 4π λ Θ i,n sin 2πf n t 1 + ϕ n , n = 1, 2, 3 where Θ i,n and ϕ n are the searching values of the amplitude and the initial phase of the nth sinusoidal function, respectively. In the process of compensating for each sinusoidal function separately, the intensity at f dc in the RD domain gradually increases. According to (32), when all of the sinusoidal functions are compensated for, the signal can be expressed as Transforming (36) into the RD domain yields where T CPI is the CPI and f dc is the Doppler center frequency. As shown by (37), after all three sinusoidal functions have been compensated for, the signal can be completely focused at f dc , i.e., the intensity at f dc is maximal. Correspondingly, thef n ,φ n , andΘ i,n of all three sinusoidal functions for the ith scatter point are matched. However, as shown in Section 2.2, because the amplitude of the sinusoidal function at different scatter points of a ship is not the same, the amplitudes of other scatter point need to be matched.
Similarly, Θ i,1 is searched and the 1st sinusoidal function of (32) is compensated for by searching for values in a range cell including the scatter point to be matched. When the intensity atf 1 is minimal, the searching value of the amplitude is the matched value ofΘ i,1 . Then, Θ i,2 and Θ i,3 are matched with the same method. Furthermore, when the amplitudes of all scatter points are matched, the matched values of Θ i,n , f n , and ϕ n for all scatter points can be obtained. When the parameters of rotational motion are non-constant in the whole CPI, the rotational motion cannot be matched accurately. Fortunately, the rotation should be constant in a sub-aperture, and the proposed method can be used in a sub-aperture.

Ship Imaging with BPA and Matched Parameters of Sinusoidal Functions
As mentioned above, the influence of rotational motion can be matched by the sum of three sinusoidal functions, and a matching method for the parameters of sinusoidal functions has been proposed. Although the influence of rotational motion can be matched, different from static targets, each scatter point of a ship moves differently due to the 3D rotational motion of the ship. Most of the existing imaging algorithms, such as (Range Doppler) RD, (Chirp Scaling) CS, and SPECtral ANalysis (SPECAN), cannot be used to image scatter points with different motions.
After range compression, the echo signal is distributed in different range cells. In the BPA, the signal is integrated along its slant range, and the integration result is the image of the target. The BPA calculates the range migration of different scatter points separately and can be used as the imaging method of a ship. As shown in (16), the imaging for the ith scatter point can be expressed as 2r where s rm is the echo model after the range compression and expressed as (16);r Ri (t) is the matched slant range of the ith scatter point and can be written aŝ

Simulation Results and Discussion
In this Section, numerical experiments are carried out to verify the effectiveness of the proposed multiple sinusoidal function matching-based rotational ship imaging method in GEO SAR. In Section 4.1, the matching method for the frequencies of sinusoidal functions, which was proposed in Section 3.2, is verified. Then, in Section 4.2, the matching method for the amplitudes and initial phases of sinusoidal functions, which was proposed in Section 3.3, is verified. Then, in Section 4.3, based on the matched values in Sections 4.1 and 4.2, the BPA for a rotating target, which was proposed in Section 3.4, is verified. Finally, in Section 4.4, the 2D raw echoes of a rotating ship are generated in the apogee of GEO SAR to verify the imaging method proposed in Section 3.1. The system parameters of GEO SAR used for the simulation are listed in Table 1; the parameters of rotational motion for the simulation are listed in Table 2.  Table 2. The parameters of the rotational motion of the ship.

Simulation Results for Matching the Frequencies of Sinusoidal Function
Based on (8), we generated the echoes of a ship target represented by a scatter point b i p i z i = −100 100 0 (m) with rotational motion into the background raw data of GEO SAR. The parameters of the rotational motion are shown as Table 2. Furthermore, the imaging of a ship just by existing BPA, which can be obtained by (38) and (39) without the matched results for rotational motion, is shown in Figure 6. From Figure 6, although the range migration and phase caused by the GEO platform has been compensated, a ship cannot be focused on without rotational motion compensation, and the signal energy has spread in a range cell.

 
SAR. The parameters of the rotational motion are shown as Table 2. Furthermore, the imaging of a ship just by existing BPA, which can be obtained by (38) and (39) without the matched results for rotational motion, is shown in Figure 6. From Figure 6, although the range migration and phase caused by the GEO platform has been compensated, a ship cannot be focused on without rotational motion compensation, and the signal energy has spread in a range cell.  After range compression, the echo signal is as shown in Figure 7(a). Due to the motion of GEO SAR, the echo signal is distributed across several range cells. Based on the method proposed in Section 3.2, after GKT and GDP, the range migration and the phase caused by the motion of the GEO platform are compensated, and the simulation results are as shown in Figure 7 As shown by (34), after the GKT and GDP, the rotational motion can be matched by the sum of three sinusoidal functions, and the signal focuses at the frequencies of three sinusoidal functions in the RD domain [52]. Therefore, an Fast Fourier transform (FFT) that was performed along the azimuth direction and the Doppler profile are shown in Figure 8. The signal is focused at the frequencies of three sinusoidal functions, and the matched values of three frequencies ˆn f are 0.100Hz , 0.056Hz , and 0.144Hz .

Simulation Results for Matching the Amplitudes and Initial Phases of the Sinusoidal Function
In Section 4.1, the matched values of three frequencies of the sinusoidal function were obtained. As shown by (34), after the GKT and GDP, the rotational motion can be matched by the sum of three sinusoidal functions, and the signal focuses at the frequencies of three sinusoidal functions in the RD domain [52]. Therefore, an Fast Fourier transform (FFT) that was performed along the azimuth direction and the Doppler profile are shown in Figure 8 After range compression, the echo signal is as shown in Figure 7(a). Due to the motion of GEO SAR, the echo signal is distributed across several range cells. Based on the method proposed in Section 3.2, after GKT and GDP, the range migration and the phase caused by the motion of the GEO platform are compensated, and the simulation results are as shown in Figure 7 As shown by (34), after the GKT and GDP, the rotational motion can be matched by the sum of three sinusoidal functions, and the signal focuses at the frequencies of three sinusoidal functions in the RD domain [52]. Therefore, an Fast Fourier transform (FFT) that was performed along the azimuth direction and the Doppler profile are shown in Figure 8. The signal is focused at the frequencies of three sinusoidal functions, and the matched values of three frequencies ˆn f are 0.100Hz , 0.056Hz , and 0.144Hz .

Simulation Results for Matching the Amplitudes and Initial Phases of the Sinusoidal Function
In Section 4.1, the matched values of three frequencies of the sinusoidal function were obtained. As shown in Figure 8, because the intensity at a frequency of 0.100Hz is maximal, the amplitude

Simulation Results for Matching the Amplitudes and Initial Phases of the Sinusoidal Function
In Section 4.1, the matched values of three frequencies of the sinusoidal function were obtained. As shown in Figure 8, because the intensity at a frequency of 0.100 Hz is maximal, the amplitude and the initial phase of the sinusoidal function with a frequency of 0.100 Hz were matched first. The phase compensation was performed using (35). When the intensity atf 1 = 0.100 Hz was minimal (shown in Figure 9a), the matched values of the amplitude and the initial phase were equal toΘ i,1 = 0.094 and ϕ 1 = 2.618, respectively. Then, similarly, the amplitude and the initial phase of the sinusoidal function were matched with a frequency off 2 = 0.056 Hz. The intensity atf 2 = 0.056 Hz was minimal, and the simulation results are shown in Figure 9b. Finally, the 3rd sinusoidal function was matched with a frequency off 3 = 0.144 Hz, and the simulation results are shown in Figure 9c. The matched values of the amplitude and the initial phase for each sinusoidal function are shown in Table 3.
Remote Sens. 2020, 12, x FOR PEER REVIEW 19 of 27 was matched with a frequency of 3 =0.144 f Hz , and the simulation results are shown in Figure 9 (c).
The matched values of the amplitude and the initial phase for each sinusoidal function are shown in Table 3.
From Figure 9, when all of three sinusoidal functions have been matched, the intensity at dc f is largest.   Table 3) are obtained. The imaging of a target was simulated using (38), and the simulation results are shown in Figure 10. The range and azimuth profile of the imaging result are shown in Figure 11. The Peak Side Lobe Ratio (PSLR) of the target in  From Figure 9, when all of three sinusoidal functions have been matched, the intensity at f dc is largest.

The Simulation Results for the BPA
In Sections 4.1 and 4.2, the parameters of the sinusoidal functions were matched. In this section, the images with the matched values (listed in Table 3) are obtained. The imaging of a target was simulated using (38), and the simulation results are shown in Figure 10. The range and azimuth profile of the imaging result are shown in Figure 11. The Peak Side Lobe Ratio (PSLR) of the target in the range and azimuth directions was -13.27 and -13.25 dB, respectively. The simulation results show that before matching the rotational motion, the target cannot be focused (shown as Figure 6). After matching the influence of rotation, as described in in Sections 4.1 and 4.2, the target can be completely focused with the proposed imaging method.

The Simulation Results for the Proposed Rotating Ship Imaging Method
In this section, the effectiveness of the imaging methods proposed in Section 3.1 is verified. A ship is represented by five scatter points and is added to an echo signal of GEO SAR. The coordinates of five scatter points are listed in Table 4; the parameters of rotation are the same as those shown in Table 2. Due to the influence of rotation, scatter points A to D cannot be focused, which is shown in Figure 12(a). scatter point E is located at the gravity center of a ship, where the rotational motion does not influence the focus.

Specifications Value
Coordinate system OBPZ , m

The Simulation Results for the Proposed Rotating Ship Imaging Method
In this section, the effectiveness of the imaging methods proposed in Section 3.1 is verified. A ship is represented by five scatter points and is added to an echo signal of GEO SAR. The coordinates of five scatter points are listed in Table 4; the parameters of rotation are the same as those shown in Table 2. Due to the influence of rotation, scatter points A to D cannot be focused, which is shown in Figure 12(a). scatter point E is located at the gravity center of a ship, where the rotational motion does not influence the focus.

Specifications Value
Coordinate system OBPZ , m

The Simulation Results for the Proposed Rotating Ship Imaging Method
In this section, the effectiveness of the imaging methods proposed in Section 3.1 is verified. A ship is represented by five scatter points and is added to an echo signal of GEO SAR. The coordinates of five scatter points are listed in Table 4; the parameters of rotation are the same as those shown in Table 2. Due to the influence of rotation, scatter points A to D cannot be focused, which is shown in Figure 12a. scatter point E is located at the gravity center of a ship, where the rotational motion does not influence the focus.   Table 5. Based on Section 3.4, the imaging results of a rotating ship are shown in Figure 12(b). The simulation results show that without matching the rotation effect, the ship cannot be focused. Correspondingly, a ship can be completely focused when the method proposed in Section 3.1 is used. In this paper, each scatter point needs to match with three amplitudes of the sinusoidal functions and the amplitudes are matched one by one. In numerical experiments, the search accuracy of each amplitude for a scatter point is 1/1000 rad and the search scope is π π −    / 2 / 2 rad-it takes several minutes for each scatter point to be matched.

Conclusion
Compared with other SAR platforms, the CPI of GEO SAR is very long (in the order of hundreds of seconds). The existing methods, which match the influence of rotational motion within seconds of the CPI via the high-order Taylor series, are obviously inapplicable.
To solve the problems for rotating ship imaging with a long CPI of GEO SAR, a novel ship imaging method using multiple sinusoidal functions to match the rotation effect was proposed. First, the analysis results show that the influence of the rotation effect can be matched using the sum of  Table 5. Based on Section 3.4, the imaging results of a rotating ship are shown in Figure 12b. The simulation results show that without matching the rotation effect, the ship cannot be focused. Correspondingly, a ship can be completely focused when the method proposed in Section 3.1 is used. In this paper, each scatter point needs to match with three amplitudes of the sinusoidal functions and the amplitudes are matched one by one. In numerical experiments, the search accuracy of each amplitude for a scatter point is 1/1000 rad and the search scope is −π/2 π/2 rad-it takes several minutes for each scatter point to be matched.

Conclusions
Compared with other SAR platforms, the CPI of GEO SAR is very long (in the order of hundreds of seconds). The existing methods, which match the influence of rotational motion within seconds of the CPI via the high-order Taylor series, are obviously inapplicable.
To solve the problems for rotating ship imaging with a long CPI of GEO SAR, a novel ship imaging method using multiple sinusoidal functions to match the rotation effect was proposed. First, the analysis results show that the influence of the rotation effect can be matched using the sum of three sinusoidal functions, and the signal model of a rotating ship in GEO SAR was obtained. Then, the multiple sinusoidal function matching-based ship imaging method was proposed. Its procedure includes (1) the modification of the GKT to correct the range migration caused by the motion of the GEO platform. After the phase caused by the platform motion is compensated for by the GDP, the signal focuses at the frequencies of sinusoidal functions in the RD domain, and the frequencies can be matched. (2) The amplitudes and initial phases of sinusoidal functions are matched via phase compensation with parameter searching. (3) Based on the matched parameters of sinusoidal functions, the rotating ship is imaged via BPA.
Finally, numerical experiments were conducted to demonstrate the effectiveness of the proposed method in GEO SAR with an orbit inclination of 53 • . The simulation results show that the proposed method can realize rotating ship imaging with a range and azimuth PSLR of about −13.27 and −13.25 dB, respectively. In general, ships have a higher radar cross-section against the surrounding sea background because of metallic materials and corner reflection structures [53]. Sometimes, the clutter from the sea surface cannot be ignored, and it extends the Doppler bandwidth and influences the proposed method. A realization of the proposed method with clutter from the sea surface will be considered in future work.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A The Effectiveness of (12) and the Residual Error with Different Numbers of Sinusoidal Functions
To verify the effectiveness of (12), the rotational effect of the ship in (8) was determined by the sum of the sinusoidal functions via the nonlinear least squares method. The rotational motion parameters are listed in Table A1, where the amplitudes and frequencies of the rotational motion are given based on the typical values of a ship [47]. The sum of sinusoidal functions 1-3 shown in (12) is used to match the rotational motion of the scatter point b i p i z i = 100 50 10 (m) in the ship's coordinate system. The matched results of the phase error caused by rotation within 100 s of the CPI are shown in Figure A1. Using (12), the phase error caused by the rotation effect and the matched results with the sum of sinusoidal functions 1-3 is calculated and shown in Figure A1a. Correspondingly, the residual phase error after matching with the sum of sinusoidal functions 1-3 is shown in Figure A1b. The simulation results show that when the number of sinusoidal functions increases, the matched curve gradually approaches the rotational motion, and the residual phase error after matching is smaller.   Figure A1 shows that as the number of sinusoidal functions increases, the residual phase error after matching is smaller. Therefore, the Root Mean Square Errors (RMSEs) of the matching results versus the sum of different numbers of sinusoidal functions are given in Figure A2. The parameters of ϕ cs and ϕ g are listed in Table A1; the parameters of rotational motion are randomly generated where ME K is the number of Monte Carlo trials; NS K is the number of samples for each of the Monte Carlo trials; and θˆ and θ are the matched and real values for each sample, respectively. Figure A2 shows that as the number of sinusoidal functions increases, the matching result for rotational motion is more accurate. Once the number of sinusoidal functions is equal to three, the residual phase error after matching is minimal.   Figure A1 shows that as the number of sinusoidal functions increases, the residual phase error after matching is smaller. Therefore, the Root Mean Square Errors (RMSEs) of the matching results versus the sum of different numbers of sinusoidal functions are given in Figure A2. The parameters of ϕ cs and ϕ g are listed in Table A1; the parameters of rotational motion are randomly generated within Θ r,p,y ∈ −15 • 15 • , f r,p,y ∈ 1/15 Hz 1/5 Hz , and ϕ r,p,y ∈ 0 • 360 • ; and the scatter point is randomly generated within b i p i z i = −100 ∼ 100 −100 ∼ 100 −10 ∼ 10 (m). The number of sinusoidal functions varies from 1 to 7, and 500 Monte Carlo trials were implemented for each number. Each of the Monte Carlo trials had a PRF of 300 Hz and a CPI of 100 s, i.e., 30000 samples. The RMSE was calculated by where K ME is the number of Monte Carlo trials; K NS is the number of samples for each of the Monte Carlo trials; andθ and θ are the matched and real values for each sample, respectively. Figure A2 shows that as the number of sinusoidal functions increases, the matching result for rotational motion is more accurate. Once the number of sinusoidal functions is equal to three, the residual phase error after matching is minimal.  Figure A1 shows that as the number of sinusoidal functions increases, the residual phase error after matching is smaller. Therefore, the Root Mean Square Errors (RMSEs) of the matching results versus the sum of different numbers of sinusoidal functions are given in Figure A2. The parameters of ϕ cs and ϕ g are listed in Table A1; the parameters of rotational motion are randomly generated where ME K is the number of Monte Carlo trials; NS K is the number of samples for each of the Monte Carlo trials; and θˆ and θ are the matched and real values for each sample, respectively. Figure A2 shows that as the number of sinusoidal functions increases, the matching result for rotational motion is more accurate. Once the number of sinusoidal functions is equal to three, the residual phase error after matching is minimal.

Appendix B The Influence of Rotation Effect on the Range Cell Migration
As shown in (15), because r i (t) is unknown, the influence of rotational motion on range migration is also unknown. The rotation of a ship involves oscillatory motion around the equilibrium point.
Furthermore, the resolution of GEO SAR (optimal ρ r_best = 5 m) is not high compared with that of the Low Earth Orbit SAR (LEO SAR). The simulation of the influence of rotational motion on the slant range is shown in Figure A3. The rotational motion parameters, except Θ r,p,y , are listed in Table A1; the parameters of the positions of the scatter points in a ship are listed in Table A2.

Appendix B：The Influence of Rotation Effect on the Range Cell Migration
As shown in (15), because ( ) i r t is unknown, the influence of rotational motion on range migration is also unknown. The rotation of a ship involves oscillatory motion around the equilibrium point. Furthermore, the resolution of GEO SAR (optimal ρ = _ 5 r best m ) is not high compared with that of the Low Earth Orbit SAR (LEO SAR). The simulation of the influence of rotational motion on the slant range is shown in Figure A3. The rotational motion parameters, except Θ , , r p y , are listed in Table  A1; the parameters of the positions of the scatter points in a ship are listed in Table A2. Table A2. The parameters of the positions of scatter points.

Specifications Value
The ship's coordinate system OBPZ , m  , i.e., the influence on slant range is less than half of the range resolution, range cell migration will not occur. Furthermore, when a scatter point is farther away from the gravity center, the influence of rotation is greater. So, when the amplitude of 3D rotation is Θ =°, , 5 r p y , °15 , and °25 , the influence on scatter point A are as shown in Figure A3(b). The simulation results show that a larger amplitude of rotation has a greater influence on the slant range. The influence on the slant range of scatter point A is less than ρ _ / 2 r best , i.e., when Θ , , r p y is less than °25 , the influence of 3D rotational motion in GEO SAR does not cause range cell migration for a scatter point whose distance from the gravity center is not further away than scatter point A. In general, it can be considered that rotational motion does not cause range cell migration in GEO SAR.  Figure A3a shows that when all amplitudes of 3D rotation are Θ r,p,y = 15 • , the influence on the slant range of scatter points A, B, and C is less than ρ r_best /2, i.e., the influence on slant range is less than half of the range resolution, range cell migration will not occur. Furthermore, when a scatter point is farther away from the gravity center, the influence of rotation is greater. So, when the amplitude of 3D rotation is Θ r,p,y = 5 • , 15 • , and 25 • , the influence on scatter point A are as shown in Figure A3b. The simulation results show that a larger amplitude of rotation has a greater influence on the slant range. The influence on the slant range of scatter point A is less than ρ r_best /2, i.e., when Θ r,p,y is less than 25 • , the influence of 3D rotational motion in GEO SAR does not cause range cell migration for a scatter point whose distance from the gravity center is not further away than scatter point A. In general, it can be considered that rotational motion does not cause range cell migration in GEO SAR.