Multiple RFI Sources Location Method Combining Two-Dimensional ESPRIT DOA Estimation and Particle Swarm Optimization for Spaceborne SAR

Abstract

Synthetic aperture radar (SAR) systems are susceptible to radio frequency interference (RFI). The existence of RFI will cause serious degradation of SAR image quality and a huge risk of target misjudgment, which makes the research on RFI suppression methods receive widespread attention. Since the location of the RFI source is one of the most vital information for achieving RFI spatial filtering, this paper presents a novel location method of multiple independent RFI sources based on direction-of-arrival (DOA) estimation and the non-convex optimization algorithm. It deploys an L-shaped multi-channel array on the SAR system to receive echo signals, and utilizes the two-dimensional estimating signal parameter via rotational invariance techniques (2D-ESPRIT) algorithm to estimate the positional relationship between the RFI source and the SAR system, ultimately combines the DOA estimation results of multiple azimuth time to calculate the geographic location of RFI sources through the particle swarm optimization (PSO) algorithm. Results on simulation experiments prove the effectiveness of the proposed method.

1. Introduction

As an active microwave remote sensing device, synthetic aperture radar (SAR) has achieved excellent performance in resource exploration, military reconnaissance, and other earth observation missions. However, the electromagnetic environment in which SAR works is becoming more complex, its working frequency bands are occupied by radio devices such as broadcasting and ground radar. For the SAR system, the signals emitted by these devices are generally called radio frequency interference (RFI) [1,2]. RFI usually causes bright lines and spots in SAR images, which dramatically degrades SAR image quality and masks targets of interest, severely weakens the performance of SAR [1,2,3,4,5]. Moreover, since strong RFI will cause a significant reduction in the signal-to-interference power ratio of the echo, we cannot accurately estimate the required radar parameters from it [2]. It is generally considered that RFI is a major challenge for SAR systems and RFI suppression becomes a necessary research topic [1,2,3,4,5,6,7].

In order to solve the RFI contamination problem, there are two main solutions to suppress RFI signals [3,4,5,6,7,8,9,10,11,12,13]. The first method is to use signal processing technology to find out the difference between RFI signal and the echo signal in the time or frequency domain, and design a filter to mitigate RFI [3,4,5,6,7]. In general, signal processing methods only utilize the signal characteristics of RFI to estimate and suppress interference. However, it is worth noting that this method can only achieve a good suppression effect in the case of accurately estimating some signal characteristics of RFI. Another method is to locate the RFI source based on the contaminated echo signal, and turn off the interference emitters or take advantage of beamforming technology to achieve spatial filtering [9,10,11,12,13]. This method mainly achieves the purpose of interference suppression by reducing the RFI signal entering the receiver. For the SAR system, taking into account the spatial characteristics of the interference signal, advanced SAR systems have the ability to adaptively adjust the antenna pattern for RFI mitigation, relying on the direction-of-arrival (DOA) of the interference signal, whose feasibility has been studied in some research [9,10]. Thus, it is greatly essential to get an accurate location of the RFI source for implementing this approach.

In the research on RFI location, the works focusing on DOA estimation were studied [7,11,12,13]. Among these research studies, References [7,11] studied the effective RFI source location methods for SAR systems. In [7], SAR images were used to record the azimuth time when the interference occurs, and on this basis, the location of RFI sources was estimated based on the scene location and the signal acquisition time. However, this method does not use the characteristics of the RFI signal and can only roughly locate the RFI. In Reference [11], a uniform linear array was designed to receive the signal, and then the multiple signal classification (MUSIC) algorithm was used for the array signal to realize a one-dimensional estimation of the direction of the RFI signal. Nevertheless, this method cannot obtain the two-dimensional angle of the RFI signal direction in the azimuth and elevation, so the two-dimensional pattern of the antenna cannot be designed based on the result to achieve RFI suppression, and the precise geographic location of the RFI source cannot be calculated. Besides, References [12,13] researched the RFI location for synthetic aperture interferometric radiometry (SAIR). Reference [12] used the MUSIC algorithm to achieve high-precision two-dimensional RFI DOA estimation. However, first of all, this work was proposed for SAIR, and the structure of the Y-shaped array designed in [12] is quite complicated for the SAR system. It is the virtual hexagonal array composed of 931 array elements, which is not suitable for SAR systems. Secondly, the two-dimensional MUSIC algorithm needs to use spatial spectrum search to complete DOA estimation, the computational complexity of this process is very high, which will greatly increase the computational burden of the system in practical applications [14,15]. Finally, this method only estimates the DOA of the RFI signal and does not calculate the geographic location of the RFI source, but the location information is important for detecting ground interference sources. The literature [13] is also a study on SAIR, the same as [12]. In addition, it introduces the spatial sparsity of RFI and considers that the spatial distribution of RFI sources obeys a prior Laplace probability function. This assumption is reasonable for the RFI source in SAIR, but there is no research showing that the RFI source in SAR also obeys similar distribution characteristics.

In order to calculate the exact geographic location of multiple RFI sources of SAR with taking into account the two-dimensional DOA estimation requirements and the complexity of the system calculation, this paper presents a new method for the location of multiple independent RFI sources combining two-dimensional estimating signal parameter via rotational invariance techniques (2D-ESPRIT) and the particle swarm optimization (PSO) algorithm. The main idea of this method is to employ array signal processing technology to estimate the spatial geometric relationship between the SAR system and the RFI sources at each azimuth time, and then use the particle swarm optimization method to solve the optimization problem in the geographic location estimation of multiple independent interference sources. This method is different from other methods in three main aspects. First of all, the proposed method makes use of an L-shaped channel array to realize the two-dimensional estimation of RFI DOA, its structure is simple and effective. The reason for choosing an L-shaped array is given in Section 2.1. Moreover, compared with the MUSIC algorithms, the two-dimensional ESPRIT used in this method does not need to construct and search the spatial spectrum of the signal, thereby reducing the computational burden and improving the computational efficiency of the DOA estimation method, which is very essential when processing large-scale SAR data. Last but not the least, the method of estimating the geographic location of RFI sources proposed in this paper can accurately indicate the specific geographic coordinates of RFI sources, instead of just getting the direction of RFI sources, which can provide location information for shutting down illegal RFI sources and help the SAR system avoid interference areas.

The rest of this paper is organized as follows. Section 2 describes the geometry model of the SAR system and RFI sources, the multi-channel array model, and the signal model. In Section 3, we present the multiple RFI sources geolocation method based on the 2D-ESPRIT and PSO algorithm under the condition of an L-shaped channel array. The simulation experiment given in Section 4 proves the feasibility of the proposed method, followed by the discussion and conclusion in Section 5 and Section 6.

2. Geometry and Signal Model

2.1. Geometry Model

For the SAR system, most RFI sources are related to electronic equipment on the surface of the earth [16]. Therefore, this paper establishes a simulation model of the normal working spaceborne SAR being simultaneously interfered with by multiple RFI sources on the ground. In the process of data acquisition, the SAR satellite moves in an elliptical orbit around the earth, and the antenna points to the imaging area on the side of the flight direction. Multiple RFI sources on the ground may exist in or outside the imaging area, which may make RFI signals enter the receiver from the main lobe of the antenna, or from the side lobes. Figure 1a shows the geometric model between the SAR system and multiple RFI sources when acquiring real echo signals and being interfered with by RFI signals in the east-north-up (ENU) coordinate system. It should be noted that, unless indicated, all geographic coordinates in this article are in the ENU coordinate system. The light yellow area on the ground shows the main beam coverage area of the SAR system antenna at the azimuth time of $t_a$. As shown in this figure, there are some RFI sources on the ground that continuously transmit interference signals to the SAR system, and when the azimuth time is $t_a$, the distance between the k-th interference source and the SAR system can be expressed as $R^k\left(t_a\right)$. During the data acquisition process, with the relative movement between the SAR system and RFI sources, the distance between them changes with the azimuth time. Figure 1b shows that the arrangement of the multi-channel array of the SAR system in Figure 1a.

As can be seen from Figure 1b, (2M-1) channels are arranged at equal distance on the Y-axis and the Z-axis of the antenna coordinate system, and there are M channels on each axis and the distance between adjacent channels is d. Among them, the channel located at the origin of the coordinate belongs to the Y-axis and Z-axis arrays at the same time, and it is used to transmit and receive signals, while other channels only receive signals. Thus, at each azimuth time, the SAR system transmits signals through channel 1, and all channels receive echo signals. Moreover, we define the elevation angle ${\theta }^k$ and the azimuth angle ${\phi }^k$ of the k-th RFI signal as shown in Figure 1b. The purpose of DOA estimation is to estimate the value of these two angles.

The reason for choosing the L-shaped multi-channel array is as follows. First of all, different from the study of single interference source localization [17], in order to estimate the two-dimensional DOA of the interference signal at each azimuth time and apply it to multiple RFI sources geolocation, deploying multiple receiving channels as a non-linear array on the SAR system is necessary. Furthermore, for the SAR system, the accurate estimation of the geographic coordinates of the RFI source must be achieved based on multiple azimuth time DOA estimation results. Therefore, this method requires the system to have a good estimation accuracy of the direction of the RFI signal at each azimuth time. In addition, if the calculation complexity required for DOA estimation at each azimuth time is too high, the calculation burden of the system will be very heavy. Among various array models, the L-shaped array not only has the advantages of simple structure and low system design difficulty, but also it can offer a larger array aperture when the number of channels and the distance between channels is constant, which makes it have a higher DOA estimation accuracy [14,18]. Moreover, when processing the signal of the L-shaped array, we can simplify the two-dimensional DOA estimation problem into two independent one-dimensional DOA problems to solve [14,15], which will reduce the computational burden. Hence, considering the requirement of DOA estimation and the advantages of the L-shaped array, this paper proposes to arrange the multiple channels of the SAR system into an L-shaped to receive signals.

Besides, in order to achieve DOA estimation, the channels in the array should meet:

• All the channels are independent and isotropic.
• The distance between adjacent channels is not more than half the wavelength of the RFI signal.

2.2. Signal Model

According to the geometric model in Figure 1a, a SAR platform travels at a speed of $V$, and continuously transmits a chirp signal and receives echo signals within the azimuth time. When the SAR system is interfered with, its echo signal is composed of the real target echo signal, RFI signal, and noise signal [1,2,3,4,5,6,7,8,9]. Thus, as the range time is $t_r$ and the azimuth time is $t_a$, and $K$ RFI sources exist in the scene, the interfered echo signal $s\left(t_r,t_a\right)$ received by each channel can be written as:

$$s\left(t_r,t_a\right)={\displaystyle \sum _{k=1}^K{j}^k\left(t_r,t_a\right)}+e\left(t_r,t_a\right)+n\left(t_r,t_a\right),$$

(1)

where $j^k\left(t_r,t_a\right)$ denotes the interference signal from the k-th RFI source, $e\left(t_r,t_a\right)$ and $n\left(t_r,t_a\right)$ denote the real target echo signal and noise signal, respectively.

Due to the complex source of RFI, it is difficult to model it with a single signal form [2,3,4]. In terms of signal model, the form of sinusoidal modulation signal and linear modulation signal is simple and often used as the basic form of microwave radiation, which can be used to model RFI [2,3,4,5,6,7,8]. For the type of RFI, researchers usually classify RFI into narrowband interference (NBI) and wideband interference (WBI) according to the ratio of the RFI bandwidth to the echo signal bandwidth [2,3,4,5,6,7,8]. Furthermore, since the source of WBI is generally other radar systems, it can be modeled by the linear modulation signal whose parameters change over time [8], which is called wideband linear frequency modulated interference (WBLFMI) in this article. Meanwhile, NBI usually comes from communication equipment or television facilities with narrow bandwidth [4,5,6]. Compared with the bandwidth of SAR, this type of RFI signal may occupy a very narrow bandwidth in the frequency domain [8], or even occupy only one frequency bin like a sinusoidal signal [2]. Therefore, this article further divides NBI into narrowband linear frequency modulated interference (NBLFMI) and radio frequency noise interference (RFNI). This modeling method conforms to the signal characteristics of RFI in actual situations [2,3,4,5,6,7,8] and has achieved good results in the research on interference detection and suppression [19,20]. In addition, the pulse width of most RFI signals is shorter than the echo window time of SAR, which results in the RFI observed in actual data being pulsed [7,21]. Taking these factors into consideration, this article models the RFI signal as the pulsed signal, and these three types of interference are written as:

$$\{\begin{array}{l}{j}_{RFNI}\left(t_r{}^{\prime },t_a\right)=A_i\left(t_r{}^{\prime },t_a\right)rect\left(\frac{t_r{}^{\prime }}{T_i}\right)\exp \left\{i\left[2\pi f_c{t}_r{}^{\prime }+2\pi f_i{t}_r{}^{\prime }+{\varphi }_i\left(t_a\right)\right]\right\}\\ j_{NBLFMI}\left(t_r{}^{\prime },t_a\right)=A_n\left(t_r{}^{\prime },t_a\right)rect\left(\frac{t_r{}^{\prime }}{T_n}\right)\exp \left\{i\left[2\pi f_c{t}_r{}^{\prime }+\pi k_n{t}_r{{}^{\prime }}^2+{\varphi }_n\left(t_a\right)\right]\right\}\\ j_{WBLFMI}\left(t_r{}^{\prime },t_a\right)=A_w\left(t_r{}^{\prime },t_a\right)rect\left(\frac{t_r{}^{\prime }}{T_w}\right)\exp \left\{i\left[2\pi f_c{t}_r{}^{\prime }+\pi k_w{t}_r{{}^{\prime }}^2+{\varphi }_w\left(t_a\right)\right]\right\}\\ t_r{}^{\prime }=t_r-\frac{R^k\left(t_a\right)}{C}\end{array}$$

(2)

where $A_i\left(·\right)$, $A_n\left(·\right)$, and $A_w\left(·\right)$ denote the amplitude of these three types of RFI signals, $T_i$, $T_n$, and $T_w$ are their pulse width,$f_c$ is the carrier frequency, $f_i$ is the frequency of RFNI,$k_n$ and $k_w$ are the chirp-rate of NBLFMI and WBLFMI, ${\varphi }_i\left(t_a\right)$, ${\varphi }_n\left(t_a\right)$, and ${\varphi }_w\left(t_a\right)$ denote the signal phase, $R^k\left(t_a\right)$ is the distance between the k-th RFI source and the SAR system at the azimuth time of $t_a$. In addition, it needs to be pointed out that the bandwidth of these three types of interference is relatively narrow compared to the bandwidth of the SAR signal [2,3,4,5].

Besides, in related research on RFI, the following assumptions are often mentioned generally:

• The real target echo signal, the RFI signal, and the noise signal are statistically independent [4,5].
• The RFI signal possesses dominant energy in the interfered echo signal, and the real target echo signal has a noise-like spectrum compared with the RFI signal, which means we can define the sum of the real target echo signal and noise signal as the equivalent additive noise [4,5,6]:

  $$g(t_r,t_a)=e(t_r,t_a)+n(t_r,t_a).$$

  (3)
• RFI signals from different interference sources are uncorrelated [9,10,11,12,13].
• The antennas of each channel of the SAR system are isotropic and not coupled to each other [9,10,11,12,13].
• The distance between RFI source and the spaceborne SAR system usually meets the far field condition [9,10,11].

According to Equations (1) and (3), the contaminated echo signal $s_{zm}$ received by the m-th channel on the Z-axis of the channel array of the SAR system can be expressed as:

$$s_{zm}(t_r,t_a)={\displaystyle \sum _{k=1}^K{j}_{zm}^k(t_r,t_a)}+g_{zm}(t_r,t_a),$$

(4)

where $j_{zm}^k(t_r,t_a)$ denotes the signal from the k-th RFI source received by this channel, $g_{zm}(t_r,t_a)$ denotes the equivalent additive noise composed of the sum of the real target echo signal and the noise signal.

At the same azimuth time, the interference signal from the k-th RFI source is received by channel $z_1$, and its received signal can be expressed as $j_1^k(t_r,t_a)$. The relationship between $j_1^k(t_r,t_a)$ and $j_{zm}^k(t_r,t_a)$ can be written as:

$$j_{zm}^k(t_r,t_a)=j_1^k(t_r-{\tau }_{zm}^k,t_a),$$

(5)

where ${\tau }_{zm}^k$ denotes the time difference between the interference signal of the k-th RFI source reaching channel $z_1$ and channel $z_m$. Since the RFI signal can be regarded as a narrowband signal compared with the echo signal, Equation (5) can be written as:

$$j_1^k(t_r-{\tau }_{zm}^k,t_a)\approx j_1^k(t_r,t_a)e^{-i2\pi f_c{\tau }_{zm}^k}.$$

(6)

Figure 2 shows a schematic diagram of the signal from the k-th RFI source being received by channel $z_1$ and channel $z_m$. Based on assumption E, the RFI signal received by these channels can be approximately regarded as a plane wave, which means that the two red lines representing the direction of the RFI signal in Figure 2 can be considered parallel.

According to the geometric relationship in Figure 2 and the definition of elevation angle in Figure 1b, the time difference ${\tau }_{zm}^k$ of the k-th RFI signal reaching the two channels can be obtained:

$${\tau }_{zm}^k=\frac{(m-1)d\cos {\theta }^k}{C},$$

(7)

where C is the light speed, ${\theta }^k$ is the elevation angle between the interference signal of the k-th RFI source and the array. Thus, combined with Equations (4) and (6), we can obtain:

$$s_{zm}(t_r,t_a)={\displaystyle \sum _{k=1}^K{j}_1^k(t_r,t_a)}{e}^{-i2\pi f_c{\tau }_{zm}^k}+g_{zm}(t_r,t_a).$$

(8)

The signal received by the M channels on the Z-axis can be expressed as a vector $S_Z$:

$$S_Z=\left[\begin{array}{l}{s}_{z1}(t_r,t_a)\\ s_{z2}(t_r,t_a)\\ ⋮\\ s_{zM}(t_r,t_a)\end{array}\right]=\left[\begin{array}{llll}{e}^{-i2\pi f_c{\tau }_{z1}^1}& e^{-i2\pi f_c{\tau }_{z1}^2}& \dots & e^{-i2\pi f_c{\tau }_{z1}^K}\\ e^{-i2\pi f_c{\tau }_{z2}^1}& e^{-i2\pi f_c{\tau }_{z2}^2}& \dots & e^{-i2\pi f_c{\tau }_{z2}^K}\\ ⋮& ⋮& \ddots & ⋮\\ e^{-i2\pi f_c{\tau }_{zM}^1}& e^{-i2\pi f_c{\tau }_{zM}^2}& \dots & e^{-i2\pi f_c{\tau }_{zM}^K}\end{array}\right]\left[\begin{array}{l}{j}_1^1(t_r,t_a)\\ j_1^2(t_r,t_a)\\ ⋮\\ j_1^K(t_r,t_a)\end{array}\right]+\left[\begin{array}{l}{g}_{z1}(t_r,t_a)\\ g_{z2}(t_r,t_a)\\ ⋮\\ g_{zM}(t_r,t_a)\end{array}\right]$$

(9)

Its compact form is:

$$S_Z=A_{Z}J+G_Z,$$

(10)

where $J={\left[\begin{array}{ccc}{j}_1^1(t_r,t_a)& \dots & j_1^K(t_r,t_a)\end{array}\right]}^T$ denotes the matrix composed of interference signals from K RFI sources received by channel $z_1$, and $G_Z={\left[\begin{array}{ccc}{g}_{z1}(t_r,t_a)& \dots & g_{zM}(t_r,t_a)\end{array}\right]}^T$ is the equivalent additive noise matrix of the M channels on the Z-axis, the steering matrix $A_Z$ is:

$$A_Z=\left[\begin{array}{cccc}{e}^{-i2\pi f_c{\tau }_{z1}^1}& e^{-i2\pi f_c{\tau }_{z1}^2}& \dots & e^{-i2\pi f_c{\tau }_{z1}^K}\\ e^{-i2\pi f_c{\tau }_{z2}^1}& e^{-i2\pi f_c{\tau }_{z2}^2}& \dots & e^{-i2\pi f_c{\tau }_{z2}^K}\\ ⋮& ⋮& \ddots & ⋮\\ e^{-i2\pi f_c{\tau }_{zM}^1}& e^{-i2\pi f_c{\tau }_{zM}^2}& \dots & e^{-i2\pi f_c{\tau }_{zM}^K}\end{array}\right].$$

(11)

In accordance with Equation (7) and the signal wavelength $\lambda =\frac{C}{f_c}$, $A_Z$ can be further transformed into:

$$A_Z=\left[\begin{array}{cccc}1& 1& \dots & 1\\ e^{-i2\pi \frac{d}{\lambda }\cos {\theta }^1}& e^{-i2\pi \frac{d}{\lambda }\cos {\theta }^2}& \dots & e^{-i2\pi \frac{d}{\lambda }\cos {\theta }^K}\\ ⋮& ⋮& \ddots & ⋮\\ e^{-i2\pi \frac{\left(M-1\right)d}{\lambda }\cos {\theta }^1}& e^{-i2\pi \frac{\left(M-1\right)d}{\lambda }\cos {\theta }^2}& \dots & e^{-i2\pi \frac{\left(M-1\right)d}{\lambda }\cos {\theta }^K}\end{array}\right].$$

(12)

Similarly, the signals received by channel $y_2$ to $y_M$ on the Y-axis can be expressed as:

$$S_Y=A_{Y}J+G_Y,$$

(13)

where $S_Y={\left[\begin{array}{ccc}{s}_{y2}(t_r,t_a)& \dots & s_{yM}(t_r,t_a)\end{array}\right]}^T$, $G_Y={\left[\begin{array}{ccc}{g}_{y2}(t_r,t_a)& \dots & g_{yM}(t_r,t_a)\end{array}\right]}^T$ is the equivalent additive noise matrix of the M channels on the Y-axis. The time difference of the k-th RFI signal to channel $y_1$ and to channel $y_m$ is ${\tau }_{ym}^k$, it can be represented as:

$${\tau }_{ym}^k=\frac{(m-1)d\sin {\theta }^k\cos {\phi }^k}{C}.$$

(14)

Similar to $A_Z$, $A_Y$ is the steering matrix of the M channels on the Y-axis, and the divergence between them is the time difference. Therefore, based on Equation (14), $A_Y$ can be written as:

$$A_Y=\left[\begin{array}{cccc}{e}^{-i2\pi \frac{d}{\lambda }\sin {\theta }^1\cos {\phi }^1}& e^{-i2\pi \frac{d}{\lambda }\sin {\theta }^2\cos {\phi }^2}& \dots & e^{-i2\pi \frac{d}{\lambda }\sin {\theta }^K\cos {\phi }^K}\\ e^{-i2\pi \frac{2d}{\lambda }\sin {\theta }^1\cos {\phi }^1}& e^{-i2\pi \frac{2d}{\lambda }\sin {\theta }^2\cos {\phi }^2}& \dots & e^{-i2\pi \frac{2d}{\lambda }\sin {\theta }^K\cos {\phi }^K}\\ ⋮& ⋮& \ddots & ⋮\\ e^{-i2\pi \frac{\left(M-1\right)d}{\lambda }\sin {\theta }^1\cos {\phi }^1}& e^{-i2\pi \frac{\left(M-1\right)d}{\lambda }\sin {\theta }^2\cos {\phi }^2}& \dots & e^{-i2\pi \frac{\left(M-1\right)d}{\lambda }\sin {\theta }^K\cos {\phi }^K}\end{array}\right].$$

(15)

Finally, the signals of all channels form the SAR echo signal matrix $S$:

$$S=\left[\begin{array}{l}{S}_Z\\ S_Y\end{array}\right]=\left[\begin{array}{l}{A}_Z\\ A_Y\end{array}\right]J+\left[\begin{array}{l}{G}_Z\\ G_Y\end{array}\right]=AJ+G,$$

(16)

where $A={\left[\begin{array}{cc}{A}_Z& A_Y\end{array}\right]}^T$ and $G={\left[\begin{array}{cc}{G}_Z& G_Y\end{array}\right]}^T$. The covariance matrix of $S$ can be written as:

$$R_S=E\left[S{S}^H\right].$$

(17)

In addition, we define the interference-to-signal ratio (ISR) as the ratio of the power of the interference signal to the sum of the power of the real echo signal and the noise signal:

$$ISR=10{\log }_{10}\left(\frac{{\displaystyle \sum {\left|J\right|}^2}}{{\displaystyle \sum {\left|G\right|}^2}}\right),$$

(18)

where ${\displaystyle \sum {\left|J\right|}^2}$ denotes the sum of RFI signal power at each azimuth time, and ${\displaystyle \sum {\left|G\right|}^2}$ denotes the sum of the real target echo signal power and the noise signal power at each azimuth time.

3. Multiple RFI Sources Geolocation Method

This section introduces the multiple RFI sources geolocation method. Firstly, the RFI detection scheme is discussed. Additionally, the minimum description length (MDL) criteria is applied to RFI source number estimation. Then, the DOA results of RFI signals at each azimuth time are calculated by performing the 2D-ESPRIT algorithm on the contaminated multi-channel SAR echo signal. Finally, on the basis of the geometric relationship between the multi-channel SAR system and RFI sources, a set of equations about the geographic location of RFI sources is established, and the results are obtained by utilizing particle swarm optimization. Figure 3 shows the process of the proposed method in summary.

3.1. RFI Detection

In our previous research, an RFI detection method for spaceborne SAR based on deep learning technology was proposed, which has achieved good results in detecting interference and identifying its type in the interfered SAR echo. More details can be found in [19,20]. With this fast and accurate method, the existence of RFI in the echo can be effectively determined in this paper. The first step of this method is to transform the echo signal into the time-frequency domain using short-time Fourier transform (STFT) to obtain a time-frequency map. Next, a trained neural network takes this time-frequency map as input and determines whether there are RFI signals in it. If RFI exists, subsequent operations of the proposed method in this paper will then be performed.

3.2. RFI Source Number Estimation

After detecting the presence of RFI signals in the echo signal, estimating the number of interference sources is a key issue for the next steps. The minimum description length (MDL) criterion is a very effective method of estimating the number of signal sources in blind source separation based on information theory [22], which applies the minimax entropy principle and does not need to set a threshold [23]. Further, the MDL criterion can estimate the number of signal sources quite accurately under the condition of a high signal-to-noise ratio [23]. Because RFI makes the ISR of the SAR echo very high, it is extremely suitable for estimating the number of RFI sources from the contaminated multi-channel SAR echo signal based on the MDL criterion. Thus, this paper applies the MDL criterion to calculate multiple sets of echo signals of the SAR system and obtains the estimated value of the number of RFI sources. Specifically, for the estimation of the number of RFI sources, the MDL criterion searches the number of RFI signals by minimizing the following equation:

$$MDL(n)=L\cdot \left(N-n\right)\cdot \ln p(n)+0.5\cdot n\cdot \left(2N-n\right)\cdot \ln L,$$

(19)

where $n$ is the number of RFI sources to be estimated, $N$ is the number of SAR channels, $L$ is the sampling number, $p(n)$ is the likelihood function, which can be written as:

$$p(n)=\frac{{\displaystyle \frac{1}{N-n}}{\displaystyle \sum _{i=n+1}^N{\lambda }_i}}{{({\displaystyle \prod _{i=n-1}^N{\lambda }_i})}^{\frac{1}{N-n}}},$$

(20)

where ${\lambda }_i$ denotes the i-th eigenvalue of the covariance matrix of SAR echo signal in Equation (17).

The procedure for estimating the number of RFI sources is to first calculate the covariance matrix $R_S$ of the SAR echo signal by using Equation (17) and compute the eigenvalues of $R_S$. Then, the likelihood function and the MDL function are constructed according to Equations (20) and (19). Finally, all the function values of the MDL function are calculated by setting $n$ in the range of $\left[1,N\right]$. The $n$ corresponding to the minimum of the function values is the estimated value of the number of RFI sources.

3.3. RFI DOA Estimation

RFI DOA estimation result can be obtained by applying the two-dimensional ESPRIT method to the array signal [24,25]. The ESPRIT algorithm obtains the incident angle information of the signal through the rotation invariance of the signal subspace of the array [24]. The L-shaped array, as shown in Figure 2, has two orthogonal uniform linear arrays. Each linear array can be regarded as composed of two sub-arrays, the channels on the Y-axis are divided into two sub-arrays: $Y_1$ and $Y_2$. Similarly, the channels on the Z-axis are divided into $Z_1$ and $Z_2$ as shown in Figure 4. The special arrangement of all the channels in the L-shaped array makes the signal received by the sub-array contain the DOA information of the RFI, which can be estimated by the two-dimensional ESPRIT method. The following is the DOA estimation process.

The Equations (10) and (13) give the signal form received by the channels on the Y-axis and the Z-axis in the L-shaped array. Further, the signals received by their sub-arrays $Y_1$, $Y_2$,$Z_1$, and $Z_2$ can be expressed as:

$$\{\begin{array}{l}{S}_{Y1}={\left[\begin{array}{lll}{s}_{y1}(t_r,t_a)& \dots & s_{y\left(M-1\right)}(t_r,t_a)\end{array}\right]}^T=A_{Y1}J+G_{Y1}\\ S_{Y2}={\left[\begin{array}{lll}{s}_{y2}(t_r,t_a)& \dots & s_{yM}(t_r,t_a)\end{array}\right]}^T=A_{Y2}J+G_{Y2}\\ S_{Z1}={\left[\begin{array}{lll}{s}_{z1}(t_r,t_a)& \dots & s_{z\left(M-1\right)}(t_r,t_a)\end{array}\right]}^T=A_{Z1}J+G_{Z1}\\ S_{Z2}={\left[\begin{array}{lll}{s}_{z2}(t_r,t_a)& \dots & s_{zM}(t_r,t_a)\end{array}\right]}^T=A_{Z2}J+G_{Z2}\end{array}$$

(21)

where $A_{Y1}$ and $A_{Y2}$ are the sub-matrices of $A_Y$, $A_{Z1}$, and $A_{Z2}$ are the sub-matrices of $A_Z$, $G_{Y1}$, $G_{Y2}$, $G_{Y2}$, and $G_{Z2}$ are the corresponding noise vectors.$A_{Y1}$ and $A_{Y2}$ have the following forms:

$$A_{Y1}=\left[\begin{array}{ccc}{e}^{-i2\pi \frac{d}{\lambda }\sin {\theta }^1\cos {\phi }^1}& \dots & e^{-i2\pi \frac{d}{\lambda }\sin {\theta }^K\cos {\phi }^K}\\ ⋮& \ddots & ⋮\\ e^{-i2\pi \frac{\left(M-2\right)d}{\lambda }\sin {\theta }^1\cos {\phi }^1}& \cdots & e^{-i2\pi \frac{\left(M-2\right)d}{\lambda }\sin {\theta }^K\cos {\phi }^K}\end{array}\right],$$

(22)

$$A_{Y2}=\left[\begin{array}{ccc}{e}^{-i2\pi \frac{2d}{\lambda }\sin {\theta }^1\cos {\phi }^1}& \dots & e^{-i2\pi \frac{2d}{\lambda }\sin {\theta }^K\cos {\phi }^K}\\ ⋮& \ddots & ⋮\\ e^{-i2\pi \frac{\left(M-1\right)d}{\lambda }\sin {\theta }^1\cos {\phi }^1}& \cdots & e^{-i2\pi \frac{\left(M-1\right)d}{\lambda }\sin {\theta }^K\cos {\phi }^K}\end{array}\right].$$

(23)

The relationship between Equations (22) and (23) can be written as:

$$\left\{\begin{array}{l}{A}_{Y2}=A_{Y1}{\Phi }_Y\\ {\Phi }_Y=diag\left(\begin{array}{ccc}{e}^{-i2\pi \frac{d}{\lambda }\sin {\theta }^1\cos {\phi }^1}& \cdots & e^{-i2\pi \frac{d}{\lambda }\sin {\theta }^K\cos {\phi }^K}\end{array}\right)\end{array}\right..$$

(24)

Similarly, the relationship between $A_{Z1}$ and $A_{Z2}$ can be presented as:

$$\left\{\begin{array}{l}{A}_{Z2}=A_{Z1}{\Phi }_Z\\ {\Phi }_Z=diag\left(\begin{array}{ccc}{e}^{-i2\pi \frac{d}{\lambda }\cos {\theta }^1}& \cdots & e^{-i2\pi \frac{d}{\lambda }\cos {\theta }^K}\end{array}\right)\end{array}\right..$$

(25)

So, according to Equations (24) and (25), Equation (21) is further derived as:

$$\left\{\begin{array}{l}{S}_{Y1}=A_{Y1}J+G_{Y1}\\ S_{Y2}=A_{Y1}{\Phi }_{Y}J+G_{Y2}\\ S_{Z1}=A_{Z1}J+G_{Z1}\\ S_{Z2}=A_{Z1}{\Phi }_{Z}J+G_{Z2}\end{array}\right..$$

(26)

It can be easily found that ${\Phi }_Y$ and ${\Phi }_Z$ are diagonal matrices consisting of angle information, and DOA can be estimated by solving them. Then, we perform pairwise covariance operations on the signals in Equation (26) to construct the covariance matrix of these signal matrices:

$$\{\begin{array}{l}{W}_1=E\left[S_{Y1}{S}_{Z1}^H\right]=A_{Y1}E\left[J{J}^H\right]A_{Z1}^H+E\left[G_{Y1}{G}_{Z1}^H\right]=A_{Y1}{R}_J{A}_{Z1}^H+R_{g1}\\ W_2=E\left[S_{Y2}{S}_{Z1}^H\right]=A_{Y1}{\Phi }_{Y}E\left[J{J}^H\right]A_{Z1}^H+E\left[G_{Y2}{G}_{Z1}^H\right]=A_{Y1}{\Phi }_Y{R}_J{A}_{Z1}^H+R_{g2}\\ W_3=E\left[S_{Y1}{S}_{Z2}^H\right]=A_{Y1}E\left[J{J}^H\right]{\Phi }_Z^H{A}_{Z1}^H+E\left[G_{Y1}{G}_{Z2}^H\right]=A_{Y1}{R}_J{\Phi }_Z^H{A}_{Z1}^H+R_{g3}\\ W_4=E\left[S_{Y2}{S}_{Z2}^H\right]=A_{Y1}{\Phi }_{Y}E\left[J{J}^H\right]{\Phi }_Z^H{A}_{Z1}^H+E\left[G_{Y2}{G}_{Z2}^H\right]=A_{Y1}{\Phi }_Y{R}_J{\Phi }_Z^H{A}_{Z1}^H+R_{g4}\end{array}$$

(27)

where $R_J=E\left[J{J}^H\right]$ is the autocorrelation matrix of RFI signal vector, $R_{g1}$, $R_{g2}$, $R_{g3}$, and $R_{g4}$ are the noise power matrices. The formulas in Equation (27) are combined into a new matrix $W$:

$$W={\left[\begin{array}{cccc}{W}_1& W_2& W_3& W_4\end{array}\right]}^T.$$

(28)

When it is considered that the signal energy of RFI is much greater than the energy of echo and noise, $R_{g1}$, $R_{g2}$, $R_{g3}$, and $R_{g4}$ in $W$ can be considered to be approximately zero. On this basis, the singular value decomposition of $W$ can get its signal subspace $D$ [14,25]:

$$D=\left[\begin{array}{c}{D}_1\\ D_2\\ D_3\\ D_4\end{array}\right]=\left[\begin{array}{c}{A}_{Y1}T\\ A_{Y1}{\Phi }_{Y}T\\ A_{Y1}{\Phi }_Z^{H}T\\ A_{Y1}{\Phi }_Y{\Phi }_Z^{H}T\end{array}\right],$$

(29)

where $D_1$, $D_2$, $D_3$, and $D_4$ are the singular value decomposition results of $W_1$, $W_2$, $W_3$, and $W_4$, respectively, and $T$ is an invertible eigenvector matrix. The relationship between them can be expressed as:

$$\left\{\begin{array}{l}{D}_3=A_{Y1}{\Phi }_Z^{H}T=\left(A_{Y1}T\right)T^{-1}{\Phi }_Z^{H}T=D_1{T}^{-1}{\Phi }_Z^{H}T=D_1\Theta \\ D_4=A_{Y1}{\Phi }_Y{\Phi }_Z^{H}T=\left(A_{Y1}{\Phi }_{Y}T\right)T^{-1}{\Phi }_Z^{H}T=D_2{T}^{-1}{\Phi }_Z^{H}T=D_2\Theta \end{array}\right.,$$

(30)

where the relationship between $\Theta$ and ${\Phi }_Z^H$ is:

$$\Theta =T^{-1}{\Phi }_Z^{H}T.$$

(31)

Further, we can solve Equation (30) and get:

$$\Theta ={\left[\begin{array}{c}{D}_1\\ D_2\end{array}\right]}^{+}\left[\begin{array}{c}{D}_3\\ D_4\end{array}\right],$$

(32)

where the symbol “+” represents the Moore–Penrose pseudo-inverse of the matrix.

It is obvious that the diagonal channels of ${\Phi }_Z^H$ are identical to the eigenvalues of $\Theta$, which can be expressed as [14]:

$${\Phi }_Z^H=T\Theta T^{-1}.$$

(33)

According to Equation (31), with the assistance of performing a singular value decomposition on $\Theta$, its eigenvalue matrix $C_Z=diag\left(\begin{array}{ccc}{C}_1& \cdots & C_K\end{array}\right)$ is the estimation of ${\Phi }_Z^H$. In accordance with Equation (25), the elevation angle ${\theta }^k$ of the k-th RFI source can be estimated:

$${\widehat{\theta }}^k=\arccos \left(\frac{angle\left(C_k\right)\cdot \lambda }{2\pi d}\right).$$

(34)

Correspondingly, according to Equations (29) and (30), we can also get the relationship between $D_1$ and $D_2$, $D_3$ and $D_4$:

$$\left\{\begin{array}{l}{D}_2=A_{Y1}{\Phi }_{Y}T=\left(A_{Y1}T\right)T^{-1}{\Phi }_{Y}T=D_1{T}^{-1}{\Phi }_{Y}T=D_1\Psi \\ D_4=A_{Y1}{\Phi }_Y{\Phi }_Z^{H}T=A_{Y1}{\Phi }_Z^H{\Phi }_{Y}T=\left(A_{Y1}{\Phi }_Z^{H}T\right)T^{-1}{\Phi }_{Y}T=D_3{T}^{-1}{\Phi }_{Y}T=D_3\Psi \end{array}\right.,$$

(35)

where ${\Phi }_Y{\Phi }_Z^H={\Phi }_Z^H{\Phi }_Y$ is used because both ${\Phi }_Y$ and ${\Phi }_Z$ are diagonal matrices.

The solution of Equation (35) is:

$$\Psi ={\left[\begin{array}{c}{D}_1\\ D_3\end{array}\right]}^{+}\left[\begin{array}{c}{D}_2\\ D_4\end{array}\right].$$

(36)

Then, the estimation of ${\Phi }_Y$ is derived from the eigenvalue matrix $E_Y=diag\left(\begin{array}{ccc}{E}_1& \cdots & E_K\end{array}\right)$ of the singular value decomposition of $\Psi$:

$${\Phi }_Y=T\Psi T^{-1}.$$

(37)

The azimuth angle ${\phi }^k$ can be estimated by combining Equations (34) and (37):

$${\hat{\phi }}^k=\arccos \left(\frac{angle\left(E_k\right)\cdot \lambda }{2\pi d\cdot \sin {\widehat{\theta }}^k}\right).$$

(38)

It is worth mentioning that Equations (33) and (37) independently perform singular value decomposition of $\Theta$ and $\Psi$, and their eigenvalues are the estimated ${\Phi }_Z^H$ and ${\Phi }_Y$. However, the order of the eigenvalues and eigenvectors obtained by the two singular value decompositions may not be correspond. It means that after getting the estimation of ${\Phi }_Z^H$ and ${\Phi }_Y$, we need to perform the pair-matching procedure before we can continue to use Equations (34) and (38) to obtain the DOA estimation. For this reason, the eigenvector matrices in Equations (33) and (37) may not be the same, they should be rewritten as:

$$\begin{array}{l}{\Phi }_Z^H=T_1\Theta T_1{}^{-1}\\ {\Phi }_Y=T_2\Psi T_2{}^{-1}\end{array}$$

(39)

where $T_1$ and $T_2$ are the eigenvectors obtained by performing singular value decomposition on $\Theta$ and $\Psi$, respectively.

In the case of a low ISR, the pair-matching relationship can be obtained by constructing the cost function of the elevation angle and the azimuth angle and using an optimization algorithm to solve it [26]. However, this method will cause a large computational burden, and for RFI, the energy of the interference in the contaminated SAR signal is much greater than the energy of the echo and noise, which means the method of using optimization algorithms is unnecessary when dealing with the RFI DOA estimation problem in the situation described in this paper. Therefore, we adopt an effective method with lower computational complexity. Firstly, we use the eigenvectors to construct the following equation, and then realize the angle pairing according to the order of the maximum value [25]:

$$Q=T_2^H{T}_1.$$

(40)

The specific approach is that after obtaining $T_1$ and $T_2$, assume that the order of $C_Y$ is correct, then calculate $Q$ based on Equation (40), and adjust the order of $E_Y$ according to the coordinates of the largest element in each row of $Q$, the pair-matched $C_Z$ and $E_Y$ are obtained in this way.

The steps of the RFI DOA estimation method used in this section are as follows.

• Step 1: In accordance with Equation (21) to get $S_{Y1}$ , $S_{Y2}$, $S_{Z1}$, and $S_{Z2}$.
• Step 2: Construct $W$ with Equation (27), and perform singular value decomposition on it to get the signal subspace $D$.
• Step 3: Calculate $\Theta$ and $\Psi$ with Equations (32) and (36).
• Step 4: Perform singular value decomposition on $\Theta$ and $\Psi$ to get $C_Z$ and $T_1$, $E_Y$ and $T_2$.
• Step 5: Perform pair-matching procedure with Equation (40).
• Step 6: Obtain ${\hat{\phi }}_k$ and ${\widehat{\theta }}_k$ with Equations (34) and (38).

3.4. RFI Sources Geolocation Estimation

By means of applying the RFI DOA estimation algorithm in Section 3.3 to the echo signal at each azimuth time when the SAR is interfered with, the positional relationship between the interference signal of K RFI sources and the radar array at each azimuth time can be calculated.

We define $U$ represents the DOA matrix at each azimuth time, and $u_i$ is DOA estimation of K RFI signals at the i-th azimuth time:

$$\left\{\begin{array}{l}U=\left[\begin{array}{cccc}{u}_1& u_2& \cdots & u_N\end{array}\right]\\ u_i=\left[\begin{array}{cccc}{u}_i^1& u_i^2& \cdots & u_i^K\end{array}\right]\\ u_i^k={\left(\begin{array}{cc}{\hat{\theta }}_i^k& {\hat{\phi }}_i^k\end{array}\right)}^T\\ 1\le i\le N\&i\in N^{*},1\le k\le K\&k\in N^{*}\end{array}\right..$$

(41)

Thus, at each azimuth time, a vector $\tilde{V}$ from the radar array to each interference source can be calculated by the DOA estimation result, its computation method approximates to Equation (41):

$$\left\{\begin{array}{l}\tilde{V}=\left[\begin{array}{cccc}{\tilde{v}}_1& {\tilde{v}}_2& \cdots & {\tilde{v}}_N\end{array}\right]\\ {\tilde{v}}_i=\left[\begin{array}{ccc}{\tilde{v}}_i^1& \begin{array}{cc}{\tilde{v}}_i^2& \cdots \end{array}& {\tilde{v}}_i^K\end{array}\right]\\ {\tilde{v}}_i^k={\left(1,\sec {\widehat{\theta }}_i^k\cdot \cot {\hat{\phi }}_i^k,\tan {\widehat{\theta }}_i^k\right)}^T\\ 1\le i\le N\&i\in N^{*},1\le k\le K\&k\in N^{*}\end{array}\right.,$$

(42)

where ${\tilde{v}}_i^k$ represents the vector direction from the SAR system to the k-th RFI source at the i-th azimuth time calculated from the RFI DOA estimation result, and $\tilde{V}$ is the matrix of these estimated vector directions of all azimuth time. The schematic diagram of solving the RFI source location is shown in Figure 5.

We define ${\Upsilon }_i={\left(x_i,y_i,z_i\right)}^T\left(1\le i\le N\&i\in N^{*}\right)$ as the spatial coordinate of the radar array at the i-th azimuth time, and $P_j^k={\left(x^k,y^k,z^k\right)}^T$ denotes the spatial coordinates of the k-th RFI source. Accordingly, the vector $V$ connecting the k-th RFI source and the radar array at all azimuth time is written as:

$$\left\{\begin{array}{l}V=\left[\begin{array}{cccc}{v}_1& v_2& \cdots & v_N\end{array}\right]\\ v_i=\left[\begin{array}{ccc}{v}_i^1& \begin{array}{cc}{v}_i^2& \cdots \end{array}& v_i^K\end{array}\right]\\ v_i^k={\Upsilon }_i-P_j^k\\ 1\le i\le N\&i\in N^{*},1\le k\le K\&k\in N^{*}\end{array}\right.,$$

(43)

where $v_i^k$ denotes the real vector direction from the SAR system to the k-th RFI source at the i-th azimuth time.

Obviously, if DOA estimates are accurate enough, $\tilde{V}$ and $V$ are parallel, which can be written as $\tilde{V}\parallel V$. At the i-th azimuth time, this parallel relationship can be expressed as:

$$\frac{x_i-x^k}{1}=\frac{y_i-y^k}{\sec {\widehat{\theta }}^k\left(i\right)\cdot \cot {\hat{\phi }}^k\left(i\right)}=\frac{z_i-z^k}{\tan {\widehat{\theta }}^k\left(i\right)}.$$

(44)

Theoretically, a straight line equation like Equation (44) can be obtained at each azimuth time, and the intersection of all straight lines is the position of the RFI source, so the position of the RFI source can be solved by combining the linear equations at all azimuth time. However, the inevitable noise and systematic errors will cause DOA estimation errors, which will make the equations difficult to solve. Hence, we propose to transform this problem into an optimization problem, which is to find a point in the space coordinate system so that the sum of the distances from the point to each straight line is the smallest. It is worth noting that since almost all RFI sources are within a certain range from the ground, we can approximate that the RFI sources exist on an ellipsoid centered at the center of the earth. This will serve as a constraint for the above optimization problem In the earth-centered, earth-fixed (ECEF) coordinate system, a point on the earth can be expressed as:

$$\left\{\begin{array}{l}{\left(P_j^k\right)}^T{\Lambda }^{-1}{P}_j^k=1\\ \Lambda =diag\left(\begin{array}{ccc}{R}_a{}^2& R_a{}^2& R_b{}^2\end{array}\right)\end{array}\right. ,$$

(45)

where $R_a$ and $R_b$ are respectively the semi-major axis and semi-minor axis of the earth, so the constraint condition is defined as:

$$\left\{\begin{array}{l}1\le {\left(H{P}_j^k\right)}^T{\Lambda }^{-1}H{P}_j^k\le {\chi }^2\\ \Lambda =diag\left(\begin{array}{ccc}{R}_a{}^2& R_a{}^2& R_b{}^2\end{array}\right)\end{array}\right. ,$$

(46)

where $H$ represents the transformation matrix from the ENU coordinate system to the ECEF coordinate system, $\chi$ represents the range ratio. In this paper, we assume that the height of the RFI source does not exceed 10 km, so it should satisfy $\chi \in \left(1,1.0016\right]$. The optimization problem for solving the coordinate position of the k-th RFI source is expressed as:

$$\left\{\begin{array}{l}\min {\displaystyle \sum _{i=1}^{N_a}{‖P_j^k-{\Upsilon }_i-\left[P_j^k-{\Upsilon }_i\right]\cdot {\tilde{v}}_i^k‖}_2}\\ s.t.\hspace{1em}1\le {\left(H{P}_j^k\right)}^T{\Lambda }^{-1}H{P}_j^k\le {\chi }^2\\ \hspace{1em}1\le i\le N\&i\in N^{*}\end{array}\right. .$$

(47)

The optimization problem in Equation (47) is non-convex and needs to be solved in its nonlinear constraints, the particle swarm optimization is applied to solve it. As a method with the ability to solve non-convex problems and global search capabilities, the PSO algorithm is easy to be implemented and the computational cost is low in memory requirements [27]. It is also a derivative-free method that can achieve convergence rapidly. Not only that, it can find a good solution in a complex search range based on the experience of particle swarms [27]. It first randomly initializes a group of particles, and then continuously searches for suitable function values in the domain. In this process, these particles share information, which enables them to constantly adjust their own optimization direction and finally find the optimal solution in the complex search space [27]. In the solution process, we define the size of the particle swarm as $N_{PSO}$ and use the optimization function and constraint conditions in Equation (47) as the fitness function and search range of the PSO algorithm. As mentioned above, the solution of this optimization problem is a three-dimensional vector, so during the t-th iteration the position of the j-th particle in the search space can be represented by a three-dimensional vector ${\rho }_j^t=\left({\rho }_{j1}^t,{\rho }_{j2}^t,{\rho }_{j3}^t\right)$, and its velocity is ${\mu }_j^t=\left({\mu }_{j1}^t,{\mu }_{j2}^t,{\mu }_{j3}^t\right)$. The shared information of the particle swarm during the t-th iteration calculation is: a) Individual optimal position: the best position that the j-th particle has ever been to is ${\delta }_j^t=\left({\delta }_{j1}^t,{\delta }_{j2}^t,{\delta }_{j3}^t\right)$; (b) Group optimal position: the global best position of the group discovered so far is expressed as ${\delta }_g^t=\left({\delta }_{g1}^t,{\delta }_{g2}^t,{\delta }_{g3}^t\right)$. The position of each particle in any iteration represents a possible solution. After each iteration is completed, fitness is used to select and update the individual optimal position and the group optimal position:

$$\left\{\begin{array}{l}{\mu }_j^{t+1}={\beta }^t{\mu }_j^t+{\upsilon }_1{\zeta }_1^t\left({\delta }_j^t-{\rho }_j^t\right)+{\upsilon }_2{\zeta }_2^t\left({\delta }_g^t-{\rho }_j^t\right)\\ {\rho }_j^{t+1}={\rho }_j^t+{\mu }_j^{t+1}\end{array}\right.,$$

(48)

where $\beta$ denotes the inertia weight, ${\upsilon }_1$ and ${\upsilon }_2$ are the cognitive and group scaling parameters and their value range is $\left(0,+\infty \right)$, ${\zeta }_1^t$ and ${\zeta }_2^t$ are two random values selected in the t-th iteration within the range $\left[0,1\right]$. According to Equations (47) and (48), the solver of the position optimization problem can be constructed.

In summary, the procedure of the RFI sources geolocation estimation method used in this section can be summarized as follow. Firstly, according to the spatial coordinate ${\Upsilon }_i$ at the i-th azimuth time of the radar array, the vector $v_i^k$ from the k-th RFI source to the radar array and the transformation matrix $H$, the optimization problem in Equation(47) can be constructed. Then, the PSO algorithm shown in Figure 6 is applied to solve it.

4. Simulation Experiments and Results

4.1. System Establishment and Parameter Setting

In Section 3, the multiple RFI sources geolocation method based on the 2D-ESPRIT and PSO algorithm has been discussed. In this section, the simulation results obtained by processing on the simulated data using the proposed method are presented. First, we establish a spaceborne multi-channel SAR system that may be affected by single or multiple RFI source. Next, on the basis of this system, the RFI transmitted from multiple different spatial locations and existing at the same time is simulated, and the contaminated echo signals received by the SAR under these conditions are obtained. The key parameters of SAR system in the simulation works are shown in

Table 1. Synthetic aperture radar (SAR) system parameters.

Parameter	Value
Pulse Width (μs)	5
Bandwidth (MHz)	90
Sampling Frequency (MHz)	100
Pulse Repetition Frequency (Hz)	1000
Receive Window Width (μs)	20
Average Orbit Height (km)	500
Average Speed (km/s)	7.61
Wavelength (m)	0.30
Distance between adjacent channel (m)	0.15
Number of channels	19
Signal-to-Noise Ratio (dB)	−15

.

For the echo signal, we take the pixel value of a real SAR image as the backscattering coefficient of the scene, in accordance with the radar and antenna parameters set in Table 1, and the relative position of SAR and the scene set as Figure 1a, the scene echo signal is calculated pulse by pulse. We have successfully applied this method in previous studies to generate echo signals [17,19,20].

For RFI sources, we define that RFI sources are independent of each other, and each RFI source is located at a place not lower than the surface of the earth, the form of interference signals emitted by them is in accordance with Equation (2). In order to fully verify the effectiveness of the proposed method, various simulation experiments are conducted. In each simulation experiment, the emission signal of the interference source is randomly generated according to Equation (2). The specific method is to first randomly generate the characteristic parameters of the signal within a certain range and then substitute these parameters into the equation to obtain the signal. In this way, the interference signal of each simulation can conform to the RFI signal characteristics and independent characteristics. The range and approach of randomly generating RFI signal parameters are shown in

Table 2. The range and approach of randomly generating RFI signal parameters.

Type of RFI	Parameters	Range	Distribution
RFNI	Frequency (KHz)	$\left[0.5,500\right]$	Uniform
	Pulse Width (μs)	$\left[1,9\right]$	Uniform
	Amplitude(dB)	$\left[-10,50\right]$	Uniform
	Phase	$\left[0,2\pi \right)$	Uniform
NBLFMI	Bandwidth (MHz)	$\left[0.09,0.9\right]$	Uniform
	Pulse Width (μs)	$\left[1,9\right]$	Uniform
	Amplitude (dB)	$\left[-10,50\right]$	Uniform
	Phase	$\left[0,2\pi \right)$	Uniform
WBLFMI	Bandwidth (MHz)	$\left[0.9,4.5\right]$	Uniform
	Pulse Width (μs)	$\left[1,9\right]$	Uniform
	Amplitude (dB)	$\left[-10,50\right]$	Uniform
	Phase	$\left[0,2\pi \right)$	Uniform

. Based on the above system parameters and interference parameters, a spaceborne SAR simulation system interfered with by RFI sources can be established.

Further, pursuant to the number of RFI sources and the type of interference signal, we design three situations when the SAR system is jammed:

• Case 1. A single RFI source.
• Case 2. Multiple RFI sources of the same interference type.
• Case 3. Multiple RFI sources of mixed types.

With the established SAR system, we simulate the SAR echo signals with different interference -to-signal ratios in these three conditions and use them as experimental data.

4.2. Simulation Experiments and Results

According to the parameters in Table 1 and Table 2, the contaminated echo signals are simulated. At each azimuth time, the interference signal of all RFI sources, the real target echo signal, and the noise signal constitute the echo signal of the SAR system according to Equation (1). Figure 7 shows the time-domain and frequency-domain forms of the demodulated SAR echo signals interfered by the three types of RFI signals.

The performance of the two-dimensional DOA algorithm is proved via the simulation experiments. Based on the obtained echo signal, we adopt the Monte Carlo simulation method for the DOA estimation method proposed in Section 3.3, and for each case, the number of Monte Carlo trials $N_{MC}$ is 1000. Generally, the root mean square error (RMSE) of the angle estimation result is defined as an index to evaluate the performance of the two-dimensional DOA algorithm:

$$RMS{E}_{DOA}=\frac{1}{K}{\displaystyle \sum _{k=1}^K\sqrt{\frac{1}{N_{MC}\cdot N_a}{\displaystyle \sum _{l=1}^{N_{MC}}{\displaystyle \sum _{i=1}^{N_a}\left[{\left({\theta }_i^k-{\widehat{\theta }}_i^k\right)}^2+{\left({\phi }_i^k-{\hat{\phi }}_i^k\right)}^2\right]}}}}.$$

(49)

In Case 1, the RFI source is located at $P_j^1={\left(0,0,973\right)}_{\mathrm{ENU}}$, and it is expressed in meters. In each experiment, we set the transmit power of the RFI source and set it to transmit one of the three interference signals, and then use the two-dimensional DOA method to process the echo signal and estimate the arrival angle. Figure 8 shows the DOA estimation error when the RFI source transmits three kinds of interference signals separately. Each point in Figure 8 represents the error between the estimated DOA and the true DOA in azimuth $\phi$ and elevation $\theta$.

For each type of RFI, we gradually increase the transmit power of the RFI source, which will make the ISR of the echo signal increase accordingly. Then, the RMSE of the DOA estimation under each ISR is calculated and recorded. The RMSE of the angle estimation for this case can be calculated by substituting the DOA estimation result into Equation (49). The RMSE of angle estimation when interference transmitting power gradually increases is shown in Figure 9.

In Case 2, three RFI sources of the same type are set in the system, and their geolocation coordinates in the ENU coordinate system respectively are $P_j^1={\left(-5.77\times {10}^4,-5.45\times {10}^4,479.82\right)}_{\mathrm{ENU}}$, $P_j^2={\left(0,0,863\right)}_{\mathrm{ENU}}$, and $P_j^3={\left(6.32\times {10}^4,5.90\times {10}^4,-388.83\right)}_{\mathrm{ENU}}$, all these coordinates are expressed in meters. As ISR is 10 dB, the angle estimation results are shown in Figure 10. The RMSE in this case is shown in Figure 11. Similar to Figure 8, each point in Figure 10 represents the error between the estimated DOA and the true DOA in azimuth $\phi$ and elevation $\theta$.

In Case 3, the number of RFI sources and geographic coordinates are the same as Case 2. As ISR is 10 dB, the angle estimation result is shown in Figure 12. The RMSE in this case and Case 2 are integrated and shown in Figure 13.

Then, we verify the performance of the PSO algorithm to calculate the spatial coordinates of the RFI source. The parameter selection refers to the principle of parameter setting in [28], $\beta$ equal to 0.5, ${\upsilon }_1$ and ${\upsilon }_2$ both equal to 2.0. Figure 14 indicates the process of using the PSO algorithm to solve the optimization problem in Equation (47).

Similarly, the RMSE of the estimated location and the true location can be used to illustrate the localization accuracy:

$$RMS{E}_{POS}=\frac{1}{K}{\displaystyle \sum _{k=1}^K\sqrt{\frac{1}{N_{MC}}{\displaystyle \sum _{l=1}^{N_{MC}}\left[{\left(x^k-{\widehat{x}}^k\right)}^2+{\left(y^k-{\widehat{y}}^k\right)}^2+{\left(z^k-{\widehat{z}}^k\right)}^2\right]}}}.$$

(50)

When the transmission power gradually increases, the RMSE of the estimated location in the three cases is shown in Figure 15.

In order to evaluate the validity of the location estimation result obtained using the proposed method, digital beamforming (DBF) is used to mitigate the RFI in the echo signal with this result as the input, and the imaging results of the suppressed signal get to illustrate the effectiveness. The RFI suppression method aligns the notch direction of the beamforming pattern with the direction of the RFI source to reduce the energy of the RFI in the echo signal [9,10]. In this paper, we use the DBF algorithm described in Reference [9]. Specifically, the estimated RFI source geographic location is utilized to calculate the elevation and azimuth angles of the RFI signal reaching the channel array at each azimuth time, and the optimal problem for RFI suppression is established according to the linearly constrained minimum variance (LCMV) method in DBF on the basis of the estimated angles [9]. The beamforming pattern weight $\omega$ is acquired by solving the optimization problem and the signal after RFI suppressing $S^{\prime }$ can be get by multiplying this weight with the channel signal $S$, which is given as [9]:

$$S^{\prime }=\omega S.$$

(51)

Combined with the results of Figure 8, Figure 10 and Figure 12, we respectively apply the estimated RFI source location results of Cases 1–3 to the DBF algorithm to obtain the beamforming pattern weights in the elevation and azimuth directions at each azimuth time. Figure 16 shows the beamforming pattern weights calculated by using the elevation and azimuth estimation results at half of the azimuth time for Cases 1–3. In Figure 16a,b, are the result of the RFI type being RFNI in Case 1, (c) and (d) are the result of the RFI type being NBLFMI in Case 2, and (e) and (f) are the result in Case 3. These figures show the weights of the beamforming pattern at various angles, and the red dashed line represents the direction of interference, and the black dashed line represents the direction of the connection between the center of the scene and the array. As shown in Equation (51), the signal in the interference direction will be weakened by multiplying these weights with the signal matrix, which achieves interference suppression.

Next, the BP imaging algorithm is applied to the signal after interference suppression to get the image results in Figure 17. Figure 17a shows the SAR image of the RFI-free echo signal. Then, the BP imaging algorithm is used to process the echo signals in Case 2 that are interfered by three RFI sources of the same type to obtain contaminated SAR images in Figure 17b,d,f. They are interfered by RFNI, NBLFMI, and WBLFMI, respectively. Figure 17c,e,g represents the imaging results after interference suppression on the signals in (b), (d), and (f). Likewise, (h) and (i) are the images for Case 3. It can be seen from the interfered images that when subjected to different types of RFI, the interference in the image is different. For the images after interference suppression, RFI is greatly suppressed, but the quality of these images is not as good as the RFI-free images. Therefore, the ISR of the signal before and after interference suppression can be used to show the degree of RFI suppression.

The ISRs before and after interference suppression are calculated after obtaining the RFI suppressed signal. The ISR before interference suppression $IS{R}_{bef}$ is obtained by Equation (18), the ISR after interference suppression $IS{R}_{aft}$ and the improvement of ISR of the signal $IS{R}_{imp}$ can be written as:

$$\begin{array}{l}IS{R}_{aft}=10{\log }_{10}\left(\frac{\left|{\displaystyle \sum {\left|S^{\prime }\right|}^2}-{\displaystyle \sum {\left|G\right|}^2}\right|}{{\displaystyle \sum {\left|G\right|}^2}}\right)\\ IS{R}_{imp}=IS{R}_{bef}-IS{R}_{aft}\end{array}$$

(52)

Table 3. The range and approach of randomly generating RFI signal parameters.

Case	Type of RFI	$\mathit{I}\mathit{S}{\mathit{R}}_{\mathit{i}\mathit{m}\mathit{p}} (\mathbf{dB})$
Case 1	RFNI	28.83
	NBLFMI	27.82
	WBLFMI	27.03
Case 2	RFNI	25.66
	NBLFMI	23.65
	WBLFMI	19.84
Case 3	Mixed RFI	21.92

shows the improvement of ISR after applying the above interference suppression methods to the signals in Figure 8, Figure 10 and Figure 12, respectively. For Case 1, because of the higher location accuracy of the RFI source, the improvement of the signal ISR is better than other situations. In Case 2, the location accuracy is the worst for the interference type WBLFMI, which also leads to the smallest improvement in its ISR. The improvement of Case 3 is better than that of WBLFMI in Case 2, which is also consistent with the result of location accuracy.

5. Discussion

In Section 4, by setting the system and interference parameters, we simulate and get the SAR echo signals contaminated by the three types of RFI, whose time domain and frequency domain forms are shown in Figure 7. Furthermore, a series of experiments are set up to discuss the performance of the 2-D ESPRIT method and the PSO algorithm under different type RFI sources and ISR conditions.

For the 2-D ESPRIT method, Figure 8, Figure 10 and Figure 12 show the DOA error between the estimation result and the true value at each azimuth time. It can be seen that the estimation error of $\theta$ does not exceed $\pm 1^{°}$, and the estimation error of $\phi$ does not exceed $\pm {0.5}^{°}$, and most of the estimation errors are small. This demonstrates that the method has a good performance in the DOA estimation of RFI sources. Comparing these three figures, we can further find that the performance of the method varies under different RFI types and ISR conditions, which can be more intuitively reflected in the RMSE results of the angle estimation. Figure 9, Figure 11 and Figure 13 show that with the gradual increase of ISR, the RMSE of the angle estimation in the three cases decrease, which illustrates that the proposed method can achieve high angle measurement accuracy when the ISR is large. It can be seen from Section 3 that the presence of noise and echo signals will make the calculation result of Equation (32) contain components that are not related to RFI, which will cause errors in the azimuth and elevation angles estimated by Equations (34) and (38). And the error increases as the ISR becomes smaller, which is consistent with the above experimental results. Meanwhile, it can be seen in the three figures that the method has different angle measurement capabilities for the three types of interference. It has the highest measurement accuracy for RFNI, followed by NBLFMI and WBLFMI. Because when the bandwidth of the RFI gradually increases, the error between the wavelength used in Equations (34) and (38) and the true wavelength of RFI increases, which will also bring error in angle measurement. In addition, we can see from the comparison between Figure 9 and Figure 11 that the optimal value of RMSE estimated for the angle of RFNI in Case 1 is ${0.0197}^{°}$, while in Case 2, the value is ${0.1523}^{°}$, which indicates that when the number of channels is the same, the DOA estimation accuracy of this method for a single RFI source is better than that of multiple RFI sources.

For the method of using DOA estimation results and PSO algorithm to estimate RFI source position coordinates, experiments also show the effectiveness of it. Illustrated in Figure 14, the particle swarm is gradually optimized to the target solution point as the iteration continues. However, the value obtained after the convergence of the PSO algorithm is slightly different from the true value, which is caused by the DOA estimation error and the random search solution process. Figure 14 also shows that the method can calculate the coordinate estimation value of the RFI source in the three cases. For Case 1, the optimal value of the RMSE of the position estimation is 284.54 m, while in Case 2 and Case 3, the value is respectively 2836.02 and 308.632 m. This proves that under the conditions of the simulation parameters set in this paper, the minimum coordinate estimation error of a single RFI source is 284.54 m, and the minimum coordinate estimation error of multiple RFI sources is 2836.02 m, so the coordinate estimation of a single RFI source is more accurate than that of multiple RFI sources. Likewise, similar to the DOA estimation result, the RMSE of the geographical location estimation in the three cases decreases with ISR increasing. Additionally, this method also has better performance for RFI signals with a narrow bandwidth. It can be proven that the location error is highly consistent with the DOA estimation error in the proposed method.

In the evaluation part, the estimated RFI source location and the DBF algorithm are utilized to suppress the RFI signal. The results in Figure 16 indicate that the obtained beamforming pattern will form a notch in a given interference direction to perform interference suppression. Although the location estimation result with errors will deviate the notch direction from the real RFI direction, the beamforming pattern also has a strong weakening effect on incoming signals within a certain error range, which means that these errors will not have a great impact on the use of the DBF algorithm to suppress interference. The results in Figure 17 and Table 3 illustrate that the RFI source location coordinates estimated by the proposed method are useful for interference suppression. They also show that the location estimation error has an impact on the interference suppression effect, but this effect is small and tolerable.

In summary, this proposed approach in this paper can estimate the geographic coordinates of the RFI source and the estimation error is acceptable.

6. Conclusions

In this paper, we have proposed a novel method for the location of multiple independent RFI sources based on the two-dimensional ESPRIT and PSO algorithm. First of all, the multi-channel spaceborne SAR echo signal model containing interference signals is established. In addition, RFI detection and interference source number estimation are performed. Furthermore, we take into account the spatial distribution characteristics of SAR signals and propose to apply the 2D-ESPRIT algorithm to the echo signal processing of the multi-channel SAR. Next, the DOA of the interference signal at each azimuth time can be achieved, which are applied to establish a set of nonlinear equations containing the coordinates of RFI sources. Finally, relying on the particle swarm optimization algorithm, we can solve these equations to obtain the geographic coordinates of RFI sources. The results obtained from the simulation experiments indicate that the proposed method provides good performance in estimating the RFI angle of arrival and geographic location. In future work, real SAR signals interfered by RFI will be used to verify the performance of this method.