A New Maritime Moving Target Detection and Tracking Method for Airborne Forward-looking Scanning Radar

Maritime moving target detection and tracking through particle filter based track-before-detect (PF-TBD) has significant practical value for airborne forward-looking scanning radar. However, villainous weather and surging of ocean waves make it extremely difficult to accurately obtain a statistical model for sea clutter. As the likelihood ratio calculation in PF-TBD is dependent on the distribution of the clutter, the performance of traditional distribution-based PF-TBD seriously declines. To resolve these difficulties, this paper proposes a new target detection and tracking method, named spectral-residual-binary-entropy-based PF-TBD (SRBE-PF-TBD), which is independent from the prior knowledge of sea clutter. In the proposed method, the likelihood ratio calculation is implemented by first extracting the spectral residual of the input image to obtain the saliency map, and then constructing likelihood ratio through a binarization processing and information entropy calculation. Simulation results show that the proposed method had superior performance of maritime moving target detection and tracking.


Introduction
The rapid development of global trade has greatly promoted the demand for marine transport and tourism, which also increases marine accidents [1][2][3].Therefore, marine search and rescue has become an important research issue and receives a lot of attention in recent years [4,5].The primary task of marine search and rescue is to discover the target under complex sea condition quickly and accurately, especially small targets, such as lifeboats carrying survivors.The airborne forward-looking scanning radar is a main detection means for marine search and rescue because of its all-weather and all-time capacity [6][7][8][9].Through the sequential scanning of radar, multi-frame data can be quickly acquired and the maritime target in the observation area can be monitored in real time.
In general, fine weather condition and calm sea surface benefit target detection and tracking because of the high signal clutter ratio (SCR) of the radar echo.In this situation, the target can be effectively detected by traditional constant false alarm rate (CFAR) technologies [10][11][12], and then the trajectory can be estimated through the tracking filter [13,14].However, villainous weather and surging of ocean waves will lead to heavy-tailed clutter and fluctuation of target scattering.In addition to the weak scattering of the target, the SCR will become quite low.In this case, traditional detection and tracking method will generate massive false alarms, making it hard to detect the true target [15].
Different from the traditional detection and tracking method, the TBD method utilizes the motion characteristics of the target and processes consecutive scans jointly to realize the fine detection and where τ is the frame sampling period, and w k represents the process noise, which obeys Gaussian distribution with mean zero and covariance Q: where q is the process noise coefficient [34].
To indicate the existence or non-existence of target in the observation area, the random variable ε k is introduced as the target existence state variable and ε k ∈ {0, 1}, where ε k =1 means there is a target in the observation area in frame k and vice versa.The above process can be modeled as a two-state Markov chain, as shown in Figure 1.In Figure 1, P birth is the probability of target birth and P death is the probability of target death, which are respectively equal to P ε k = 1|ε k−1 = 0 and P ε k = 0|ε k−1 = 1 .Thus, the two state Markov chain can be defined as

Airborne Forward-Looking Scanning Radar Measurement Model
The geometry illustration for airborne forward-looking scanning radar is shown in Figure 2.For the common overlapping area of multi-frame images, the observation area Ω k of frame k can be represented as where is the scattering coefficient of the point (x m , y n , 0) in the observation area, and N x and N y are the numbers of resolution cells along x-axis and y-axis, respectively.Assume the airplane is moving along y-axis at speed v, the waveform of transmitted signal is s (η), where η is the fast time, the scanning speed of radar beam is ω and the antenna pattern is w a (t).At slow time t, the coordinate of radar is set as (0, vt, H), and then the echo signal of the observation area in forward-looking scanning radar can be expressed as where shows the instantaneous slant range of the point (x m , y n , 0), while R x m ,y n is the initial slant range, θ x m ,y n is the initial azimuth angle and ϕ is the grazing angle.In practice, due to the narrow azimuth beam, fast scanning speed and small grazing angle [7], the slant range can be approximated as For the convenience of spatial data analysis, the time variable t and η can be represented by the spatial variables R and θ based on t = (θ − θ 0 ) ω and η = 2R/c, where θ 0 is initial scanning angle of the antenna.After the processing of matched filtering and range immigration correcting along the range dimension, the echo can be modeled as a convolution of the scattering coefficients and the two-dimensional system function [6], which is where is the two-dimensional system function, which can be written as where B is the bandwidth of the transmitted signal.
In azimuth, the resolution of radar image can be significantly enhanced by super-resolution processing, contributing to an improved two-dimensional system function [35].
After coordinate transferring of Equation ( 9), the final image in spatial domain can be represented as where h (x m , y n ) is the two-dimensional system function in spatial domain.The imaging procedure of the forward-looking scanning radar is given in Figure 2.
Assuming the size of resolution cell in the observation area is ∆ x × ∆ y , the measurement data in frame k can be expressed as The hypothesis of the existence or non-existence of a target in resolution cell (m, n) can be defined as: where I k is the amplitude of target, φ k is the random phase of target, and ς k (m, n) is the amplitude distribution of sea clutter, which is related with the sea condition and parameters of the surroundings.

The Proposed Method
In this section, the proposed SRBE-PF-TBD method is specified.The PF-TBD used in the proposed method is based on Bayesian recursion, which is derived in detail in [21].In the PF-TBD, the particles are divided into two parts: the continuing particles, which represent the continuing presence of the target, and the birth particles, which represent the newly appearing target.By using the likelihood ratio to update the mixing weights of the two part particles, these particles can be mixed into a large set, and then the target state is estimated by the large set of particles.
Therefore, the likelihood ratio calculation plays an import role in the PF-TBD.As described above, the traditional distribution-based likelihood ratio calculation is infeasible in the situation of non-stationary sea surface.The proposed method implements likelihood ratio calculation from the perspective of image information, which is independent from the prior knowledge of clutter.Based on information theory, the information of an image can be decomposed into two parts: the prior knowledge and the innovation [33].The prior knowledge part represents the statistical invariant properties of the environment, which can be considered as the redundant information and should be suppressed for target detection.The innovation part is the novel part, which contains information about target.Inspired by this rationale, the proposed method uses spectral residual calculation to remove the redundant information, and then the saliency map can be obtained.Finally, the likelihood ratio expression is constructed through a binary processing and information entropy calculation, which is named as SRBE likelihood ratio.
Based on the above, this section first presents the calculation of spectral residual binary entropy (SRBE) likelihood ratio, and then shows the implementation of the proposed SRBE-PF-TBD method.

SRBE Likelihood Ratio Calculation
Given the airborne forward-looking scanning radar image The Fourier spectrum of the image can be calculated by where F [•] means the operation of two-dimensional Fourier transform, and the corresponding amplitude spectrum and phase spectrum can be denoted as where A [•] means getting the amplitude of the input, while P [•] obtains the phase.Then, the log spectrum of the image amplitude can be expressed as Thus, the spectral residual can be obtained by where χ r ( f m , f n ) is a local mean filter, and defined by a r × r matrix: By performing two-dimensional inverse Fourier transform and Gaussian filtering, the saliency map in spatial domain can be expressed as where F −1 [•] means the operation of two-dimensional inverse Fourier transform.The Gaussian filer g (m, n) is defined by where γ is the filter parameter.Then, the normalized saliency map can be obtained by By a given binarization threshold T B , the normalized saliency map is binarized and the corresponding binary image B k is generated, which indicates the potential target area.Then, the information entropy of B k is calculated by where p 0 and p 1 , respectively, indicate the probabilities of the pixel value 0 and 1 in the binary image B k .The final SRBE likelihood ratio can be written as For the target with state x k , the coordinate of the corresponding cell can be calculated by where • denotes the operation of rounding down.The SRBE likelihood ratio corresponding to x k can be calculated by The flowchart of SRBE likelihood ratio calculation is shown in Figure 3.As shown in Figure 3, the proposed SRBE method has two key operations: one is obtaining the saliency map of the radar image by the spectral residual method, and the other is calculating the binary entropy of the normalized saliency map.Because the maritime target is generally composed of metal material and has dihedral and trihedral structures, the scattering of the target is normally higher than the scattering of the sea clutter, meaning that the target is more salient than the sea clutter in the radar image.Therefore, the saliency map can effectively indicate the potential target regions.Through the information theory, the binary entropy of the saliency map can be used to measure the cluster-degree of the saliency map.The larger the entropy value is, the smaller the cluster-degree is, indicating that the saliency map is more likely caused by the sea clutter scattering.Contrarily, the lower entropy value represents the larger cluster-degree, which means the higher probability of target existence.Therefore, the binary entropy of the saliency map can effectively indicate the confidence of target existence.Based on the above analysis, the SRBE likelihood is constructed as Equation ( 27) to update the particles, which utilizes the information of the potential target regions and the confidence of target existence in the radar image.

Implementation of SRBE-PF-TBD
By incorporating the SRBE likelihood ratio into the PF-TBD, the specific implementation of the proposed SRBE-PF-TBD is described as follows.

Calculation of birth density
A set of N b birth particles is drawn from the proposal density, whose positions are uniformly distributed within some highest measurements [22].The sampled birth particles for frame k can be represented as where q (•) is the proposal density, and the symbol "∼" means the operation of drawing samples from the proposal density.
Based on the likelihood ratio proposed in Equation ( 27), the non-normalized weights of these birth particles can be calculated by where the prior probability density Then, the normalized weights can be obtained by

Calculation of continuing density
Since the auxiliary particle filter (APF) utilizes the information of measurement in the current frame, the knowledge of current observation is incorporated into the proposal mechanism, and particles are not moved blindly in the state space, making the current target state estimated with high reliability [36].Considering these benefits, this paper adopts APF for continuing particle filtering.
A set of N c continuing particles is drawn using the dynamic equation described in Equation ( 1), which can be expressed as where p (•) represents the transition density from frame k − 1 to frame k.
The non-normalized weights of these particles are calculated with the current measurements as below.
The normalized weights of continuing particles can be calculated by bk By using weights bk (c)i , the particles , where i j N c j=1 denotes the index of particle Then, particles for frame k by Equation ( 1) and the resampled continuing particles in frame k − 1 are again drawn, which can be expressed as The non-normalized weights of newly resampled particles can be calculated by The normalized weights can be obtained by

Calculation of mixing probabilities
The mixing probability of birth density is calculated by Similarly, the mixing probability of continuing density is calculated by Then, the normalized mixing probabilities can be expressed by 4. Calculation of scaled particle weights The particle weights are scaled according to the mixing probabilities, and the scaled particle weights of birth density and continuing density are obtained: Then, the particles of birth density and continuing density are combined into one large particle set: The particles of the large set are resampled to obtain N c particles with uniform weights:

Calculation of probability of existence
The probability of existence in frame k can be expressed as

Estimation of target state
The detection threshold of existence is assumed as P th .When p ε k = 1|z 1:k > P th , the hypothesis of target existence is accepted and the target state will be estimated by

Experimental Results
To verify the proposed method, simulation experiments of maritime moving target detection and tracking based on simulated data were conducted.
The simulated data of airborne forward-looking scanning radar consisted of 30 frames, whose frame sampling period was τ = 1s, and per frame contained N x × N y = 60 × 60 bins.The target appeared in frame K b = 6, disappeared in frame K d = 21, and moved according to the dynamic equation (Equation ( 1)), while the process noise coefficient was q = 0.001.The initial state of the target was set as x 0 = 12.2 1.2 15.1 0.8 T .There were 4000 particles, where the numbers used for birth density and continuing density were N b = 2000 and N c = 2000, respectively.The initial particle velocities were sampled uniformly from interval [0, 2], and the particle positions were uniformly distributed within the 200 highest measurements.The probabilities of birth and death are defined as P birth = P death = 0.05, while the detection threshold was set as P th = 0.5.The binarization threshold T B was set to 0.5.All simulation results were averaged over 100 Monte Carlo trials.

Simulated Data
In practice, since airborne forward-looking scanning radar is used for long-range detection with small grazing angle, it is reasonable to model the sea clutter amplitude as K-distribution [29], which can be expressed as where α is the shape parameter of the distribution, while β represents the scale parameter, and is the modified second-kind Bessel function of order α − 1.
Due to the movement of the ocean, the RCS of the maritime moving target varies between radar scans.Therefore, Swerling type I model was used to model the maritime moving target in our simulation [37].The amplitude of the target can be written as where σ denotes the mean squared target amplitude.Thus, the SCR can be defined as SCR = 10log 10 σ P clutter (49) where P clutter means the mean power of the clutter in frame k and can be expressed by

Feasibility Analysis
The feasibility of the proposed method was demonstrated when SCR was set to 10dB.The true target trajectory and the estimated trajectory are shown in Figure 6.In the first two frames where the target existed, the target positions were not well estimated.However, with the accumulation of frames, the SRBE-PF-TBD method tracked the target effectively, which verified the effectiveness of SRBE-PF-TBD method.

Performance Evaluation
To analyze the performance of the proposed method, this section first gives the definition of the evaluation criteria [28,38,39], and then prsents the performance evaluation of the proposed method.

Root Mean Squared Error (RMSE) of Estimated Position
The RMSE is calculated by where M is the number of Monte Carlo experiments, and xk i , ŷk i and x k i , y k i denote the estimated position and true position of the target in the ith Monte Carlo experiment, respectively.

Average RMSE
The average RMSE is calculated by where where K num p ε k = 1|z 1:k > P th is the number of frames that satisfy the decision p ε k = 1|z 1:k > P th in the ith Monte Carlo experiment.4. False-alarm Probability

Earliest Effective Detection
Frame where F i represents the frame number that the target is first detected in the ith Monte Carlo experiment.
Figure 7 shows the probability density functions of the K-distribution corresponding to five different sets of parameters, while Figure 8 demonstrates the performance of the proposed method under these different K-distribution parameters.Figure 8a-d, respectively, shows the detection probabilities, false-alarm probabilities, earliest effective detection frame and average RMSE of the proposed method with a fixed detection threshold under different SCRs.The results in Figure 8a,c,d show that, with the increase of SCR, the detection probability increased, while the earliest effective detection frame and the average RMSE declined.In the case of the same SCR, the proposed method performed better with larger shape parameters.As for the false-alarm probabilities in Figure 8b, the false-alarm probabilities basically remained consistent under different SCRs when the distribution parameters of the K-distribution were fixed.Meanwhile, the false alarm probability decreased when the shape parameter α became larger.These are reasonable because, when the shape parameter of the K-distribution increased, the heavy-tailed effect of K-distribution became weaker.One extreme case was that, when the shape parameter tended to infinity, the K-distribution degenerated into a Rayleigh distribution.Then, the proposed method was compared with other two methods: the K-based method [28] and the ESIR method [21].Figures 9-11, respectively, show the probabilities of existence and RMSEs calculated by the three methods in cases of SCR = 8 dB, 12 dB and 16 dB.The shape parameter of K-distribution was set as α = 2 and the scale parameter aws set as β = 2 √ 2, thus P clutter = 1.In the results of probability of existence, the frames of true target existence are represented by the black dash line.From the results of probabilities of existence, the K-based method and the proposed method had better target existence estimation with higher SCR.In the case of the same SCR, the proposed method performed better than the K-based method.The ESIR method had the best estimation of target existence when there was a target.However, in the case of target non-existence, the ESIR method still gave a high probability of existence, representing a large false alarm.This is because the ESIR method is based on Rayleigh distribution, which will generate large false alarm when the clutter has the heavy-tailed feature.
In the results of RMSEs, the RMSE bound is illustrated by the black dash line.The results show that the RMSEs of these methods decreased with the increase of SCR, which means the estimated positions of the three methods gradually converged to the true position of the target with the accumulation of target existence time.In comparison, the estimation accuracy of the proposed method was superior to the other two methods.Figure 12 shows the detection probabilities of the three methods under different SCRs in the case the false-alarm probability was equal to 10 −2 .In Figure 12a, the K-distribution parameters were set as α = 2 and β = 2 √ 2. In this situation, the K-based method performed best when the SCR was below 8 dB, while the proposed method became superior to the other two methods with the increasing SCR.Meanwhile, it is worth noting that the detection probabilities of the three method were all smaller than 0.1 when the SCR was lower than 8 dB, indicating the inability to achieve reliable target detection in practical applications.In Figure 12b, the K-distribution parameters were α = 10 and β = 2 √ 2, which means a smoother K-distributed clutter compared with Figure 12a.While the SCR was lower than 5 dB in Figure 12b, the detection probability of the ESIR method was slightly higher than the other two methods, which was attributed to the approximate Rayleigh distributed clutter in this case.When the SCR was beyond 5 dB, the proposed method had better detection performance than the other methods.Based on these results, it can be concluded that the proposed method can achieve better detection performance under more extensive clutter distribution parameters compared with the other two methods.

Computational Efficiency
Computational efficiency is a key indicator for real-time processing.To compare the computational efficiency of the three methods, all experiments were done in MATLAB R2018a, while the computer CPU was Intel (R) Core (TM) i5-7500 3.40 GHz, with 16 GB RAM.
Each frame contained 60 × 60 bins, the number of Monte Carlo experiments was 300, and the average running time of single frame was calculated by where t f is the running time of single frame.The computational consumption results are given in Table 1, which shows the average running time of the three methods with different numbers of particles.The results demonstrate that the proposed method and the ESIR method were obviously superior to the K-based method in term of running time.The K-based method calculated the likelihood ratio by a numerical solution of integral operation, causing a huge computational consumption.The computational efficiency of ESIR method was slightly higher than the proposed method, but, as shown in previous experimental results, the ESIR method only worked under the Rayleigh distribution clutter and caused high false alarm under the non-Gaussian clutter background.All these experimental results demonstrate that the proposed method had fine detection and tracking performance and high computational efficiency.

Conclusions
In this paper, an effective maritime moving target detection and tracking method motivated by the human vision system is presented for airborne forward-looking scanning radar.The key idea of the proposed method is to use the saliency information of radar image to construct the SRBE likelihood ratio, which is independent from the prior knowledge of sea clutter.The proposed method effectively deals with the difficulty in target detection and tracking caused by non-Gaussian and non-stationary sea clutter.The performance of the proposed method was verified by Monte Carlo experiments using K-distribution clutter and Swerling type I target model.In the simulation experiments, the proposed method was evaluated under different K-distribution clutter parameters and compared with other two methods.The experimental results demonstrate that the proposed method had superior performance of target detection and tracking.Besides, the proposed method was computationally efficient and suitable for real-time application.
There are several continuing areas of research for improving the proposed method.The proposed method has been verified by simulated data, and future research will be to apply the proposed method to measured data.In addition, the extension of single-target situation to multi-target situation is the field to be researched.

Figure 2 .
Figure 2. Geometry illustration for forward-looking scanning radar and the flowchart of imaging processing.
Figure 4   shows the radar images of Frames 1, 5, 10, 15, 20 and 25.The target exists in Frames 10, 15 and 20, and the positions of the target are marked by red rectangles in corresponding images.As shown in Figure4, the target was difficult to distinguish from the clutter background.

Figure 4 .
Figure 4. Sampled images at different frames.

Figure 5
Figure5illustrates the distribution of particles calculated by the proposed SRBE-PF-TBD method, where the red triangles show the true positions of the target.When there was no target in Frames 1 and 5, particles were distributed according to the saliency map of the corresponding frame.In Frames 10, 15 and 20, where the target existed, the particles clustered around the target.When the target did not exist in Frame 25, the particles were distributed according to the saliency map again.

Figure 5 .
Figure 5. Distribution of particles corresponding with frames in Figure 4.

Figure 6 .
Figure 6.The estimation results of the proposed method.

Figure 7 .
Figure 7. Probability density functions of K-distribution with different shape and scale parameters.

Figure 9 .
Figure 9. Probabilities of existence and Results of RMSE in case of SCR = 8 dB.

Figure 10 .
Figure 10.Probabilities of existence and Results of RMSE in case of SCR = 12 dB.

Figure 11 .
Figure 11.Probabilities of existence and Results of RMSE in case of SCR = 16 dB.

2 √ 2 Figure 12 .
Figure 12.Detection probabilities of the three methods in case of P f a = 10 −2 .
Detection probabilities under different parameters