An Improved Adaptive Median Filtering Algorithm for Radar Image Co-Channel Interference Suppression

In order to increase the accuracy of ocean monitoring, this paper proposes an improved adaptive median filtering algorithm based on the tangential interference ratio to better suppress marine radar co-channel interference. To solve the problem that co-channel interference reduces the accuracy of radar images’ parameter extraction, this paper constructs a tangential interference ratio model based on the improved Laplace operator, which is used to describe the ratio of co-channel interference along the antenna rotation direction in the original radar image. Based on the idea of between-class variance, the tangential interference ratio threshold is selected to divide co-channel interference into high-ratio regions and low ones. Moreover, an improved adaptive median filter is used to process regions of high ratio based on the median of sub-windows, while that of low-ratio regions is processed by the adaptive median filter based on the median of current windows. Radar-measured data from Bohai Bay, China are used for algorithm validation and experimental results show that the proposed filtering algorithm performs better than the adaptive median filtering algorithm.


Introduction
Relying on the advantages of detection in low-visibility conditions, marine radar is regarded as one of the key instruments in the future of ocean monitoring, pilotage and collision avoidance [1]. Especially in terms of ocean information extraction, if reliable winds [2], waves [3], currents [4], tides [5], sea ice [6] or oil spill [7] measurements can be obtained from marine radars, the costs associated with traditional field means such as anemometers, buoys or manual visual inspection can be significantly reduced [8]. Therefore, marine radar is considered to be one of the key means of ocean information extraction today, and the land-based or shipborne marine radar have received extensive attention. However, the co-channel interference signal is of high strength, which reduces the quality of radar images and affects the accuracy of echo information extraction from the radar images. Since it is spread throughout the entire radar image, accurate suppression of the co-channel interference has become a challenge under the trend of an increasing demand for ocean monitoring.
In the processing of co-channel interference, domestic and foreign scholars have proposed many methods. Generally speaking, the methods of co-channel interference removal can be divided into hardware and software methods. The hardware approach suppresses the interference signal by adjusting the transceiver design parameters [9][10][11][12][13]. The authors of [9] proposed a coherent front-end based on the phase-synchronous optoelectronic oscillator (OEO) with low phase noise; the SNR of the echo pulses has an improvement relative to the signal before pulse accumulation. In [10], a new radio frequency interference (RFI) mitigation method based on joint fractional Fourier transform (FRFT) and complex empirical mode decomposition (CEMD) is proposed. A standard component within the radio frequency (RF) radar receiver is proposed in [11] to suppress the co-channel interference.
In [12], an improved histogram method for pulse repetition interval (PRI) deinterleaving based on pulse correlation is proposed, with the aim of radar pulse deinterleaving. However, the frequency band selection of radars is subject to many factors. The effect of the methods above is greatly reduced when all radar frequency bands are occupied. In addition, when two ships are close to each other, the effect on co-channel interference suppression is limited [13].
The software methods mainly include mean filtering, median filtering, adaptive filtering, Wiener filtering and filtering based on wavelet transform [14][15][16][17][18][19][20]. In [14], an improved adaptive mean filtering algorithm is proposed by assigning a certain weight to the pixel gray values of the noise points in the filter template, which can effectively remove pulse noise points. The authors of [15] proposed a non-local mean (NLM) filtering algorithm to enhance image quality in a highly turbid environment. In [16], a fully 3D non-local mean parallel approach is designed and implemented. To demonstrate its high applicability and scalability, different algorithm mapping strategies on a GPU architecture and multi-GPU framework are adopted. In [17], the computational model of weighted mean filtering and the characteristics of a high-performance computer architecture are studied. Moreover, an efficient hierarchical image-weighted mean filtering parallel algorithm for Open Computing Language (OpenCL) is designed and implemented. According to the respective correlation matrices of images and noise, many works have been conducted on the Wiener filter [18][19][20]. The purpose of [18] is to confirm that image quality can be improved by the median modified Wiener filter (MMWF) technique, and the authors of [19] proposed a synthetic aperture radar (SAR) image anti-speckle filter based on an extended adaptive Wiener filter (EAWF), extended guided filter (EGF) and weighted least squares (WLS) filter. In [20], contrast-limited adaptive histogram equalization (CLAHE) is used to improve contrast, and the Wiener filter is used for noise reduction. Adaptive filtering uses the filter parameter results obtained at the previous moment to automatically adjust those at the present moment to achieve optimal filtering [21][22][23][24]. In [21], the formalism for combining the Bayesian maximum entropy strategy with the variational analysis (VA) paradigm is presented to improve the Bayesian maximum entropy-variational analysis (BMEVA) performance for high-resolution radar imaging and denoising. The feasibility of integrating an adaptive filter approach for the compensation of platform motion artefacts is investigated in [22] for the extraction of respiratory motion signatures. In [23], an algorithm based on the fuzzy impulse detection technique is proposed, which can remove impulse noise efficiently from highly corrupted images while preserving image details. The authors of [24] proposed a novel adaptive Type-2 fuzzy filter for removing salt and pepper noise from images. In this study, two approaches have been proposed for finding the threshold between different types of pixels by designing a primary membership function (MF). The above software-based methods eliminate image interference or noise to a certain extent, but they still have room for improvement in image target protection or computational speed.
In this paper, we consider an improved adaptive median filtering algorithm based on the tangential interference ratio for radar co-channel interference suppression. In this work, we aim to automatically identify regions of different interference ratios to optimize the selection of filtering windows. The major contributions of this paper are as follows: • Based on the source of radar co-channel interference and the storage form of radar echo data, the Laplace operator is improved to better identify radial interference. • Radar images are processed with the improved Laplace operator and binarized to establish the tangential interference ratio model. Based on the idea of between-class variance, the tangential interference ratio threshold is determined, which provides a classification basis for different filtering methods in high-ratio regions and low ones. • We explore the improved adaptive median filtering algorithm, using different filtering windows for high-ratio regions and low ones. On the basis of protecting the details of radar echo images to the maximum extent, the gray values of interference points are replaced with the median of pixel points in the adaptive window.
The remainder of this paper is organized as follows. Section 2 introduces the improved method of Laplace operator based on the radar co-channel interference source and radar echo data storage form. Section 3 explains the principles of the OTSU algorithm firstly, and then introduces the construction process of the tangential interference ratio model and finally gives the algorithm flow. The effectiveness of the proposed algorithm is provided in Section 4, and the conclusion is given in Section 5.

Characteristics of Co-Channel Interference
During the working process of a marine radar, it will be interfered with by the electromagnetic waves emitted by other nearby radars in the same frequency band, resulting in the generation of radial interference on the radar grayscale images-that is, co-channel interference. Co-channel interference can be divided into co-channel synchronous interference and co-channel asynchronous interference.
where PRF 1 and PRF 2 denote the pulse repetition frequency transmitted by the two radars, respectively. T 1 and T 2 represent the pulse widths of the two radars, respectively; let ∆t be given by (1). When ∆t ≤ T 1 , the interference of the first radar with the other one is known as co-channel synchronous interference, and vice versa as co-channel asynchronous interference [25]. This paper mainly discusses the suppression of the latter. Considering that the radar co-channel interference signal appears as a radial bright spark with obvious directionality on the echo images, this paper improves the Laplace operator, which is suitable for point detection, to become a line detection template.

Improvement of Laplace Operator
The Laplace operator of a two-dimensional function f (x, y) is the second-order derivative defined as follows [26]: where ∇ represents the direction in which the value of a continuously differentiable multivariate function increases the fastest at a certain point. ∂ is the partial derivative of f (x, y) with respect to a variable. The difference method is used to calculate secondorder derivatives of the Laplace operator in the x and y directions, and the second-order difference in the two-dimensional function f (x, y) in the x and y directions can be obtained as follows [27]: Therefore, the difference form of the Laplace operator is It is written in the form of a matrix as follows [28]: The matrix has the same values in the four directions of up, down, left and right. Since the co-channel interferences in the situation involved in this paper are all radial bright sparks along the radial direction, the traditional Laplace operator is not suitable because of its directionality-free characteristic.
In order to better detect co-channel interference and enhance the radial grayscale mutation, the improved Laplace operator is defined as follows: The template of the improved operator can be expressed as

Improved Algorithm
The median filtering method is a nonlinear smoothing technique, which replaces the gray value of each pixel with the median value of all pixel gray values in a certain neighborhood window of the point. The principle of the median filtering method is to use a sliding template of a certain structure to sort the pixels in the template according to the size of the pixel value [29], and generate a data sequence that monotonically increases (or decreases). It is commonly used because it can remove nearly all interference or noise, with little impact on the original image. However, in addition to the large amount of computation required by the median filtering method, the original median filtering method does not have a proper distinction between interference ratios, so that it is not effective in removing co-channel interference in high-ratio regions. In this paper, the concept of the tangential interference ratio is introduced into the median filtering algorithm. The idea of between-class variance is introduced to find the threshold between high-ratio regions and low ones in the binarized image. Finally, a tangential interference ratio model is constructed to filter out radar co-channel interference.

Radar Image Collection and Analysis
The data in this paper were collected from the Bayuquan radar station in Yingkou, China. Figure 1 is an echo image received by the radar antenna. The irregular echoes on the left side of Figure 1 are the sea ice echoes, and the red box in the middle is the breakwater built to prevent the invasion of waves. As the research object of this paper, the radar cochannel interference signal appears as bright sparks throughout the direction perpendicular to the X axis, i.e., the direction of the detection radial. The two red boxes on the right of Figure 1 are examples of co-channel interference signals. The positive direction of the X axis represents the rotation direction of the antenna and the direction of the Y axis represents the direction of the detection radial. The whole image contains 1000 × 3000 pixels, the angular resolution of the X axis direction is 0.1 • and the distance resolution of the Y axis direction is 3.75 m.

Median Filtering Algorithm
The echo signal received by radar antennae is often polluted by a lot of interference and noise during its formation and transmission. In addition to the Gaussian noise generated by the machine itself, there is also co-channel interference. In order to suppress and eliminate these interferences and noise, thus improving the quality of images, it is necessary to perform denoising processing on the image-that is, filtering processing. The median filtering algorithm proposed by Tukey has been widely used to remove polluted points in images. The standard median filtering algorithm selects pixels in digital images or sequences and pixels around the adjacent pixels, and then takes the pixel value in the middle position after sorting as the pixel value of the current pixel [30]. This method allows the surrounding pixel values to be close to the true value, thereby eliminating isolated polluted points [31]. However, as the interference or noise ratio increases, the standard median filtering algorithm is not effective in preserving image details. The standard median filtering algorithm is described as follows Algorithm 1.

Algorithm 1 Median filtering algorithm
Step 1: Obtain first address of the original image and size of the image; Step 2: Open up a memory buffer to temporarily store the processing results and initialize it to 0; Step 3: Scan pixels in the binarized image through the loop statement, and sort the pixel values of each element in its neighborhood in ascending order. Finally, assign the obtained intermediate value to the pixel corresponding to the current point in the target image; Step 4: Repeat step 3 until all pixels in the original image are processed; Step 5: Copy running result from the memory buffer to the data region of the original image to complete filtering process.
In recent years, a variety of adaptive algorithms based on the median filtering algorithm have appeared. The extreme median filter (EM) proposed in [32] firstly divides all pixels into two categories according to the criterion: noise points and signal points. Then, the noise point is replaced by the median value of the neighborhood of that point according to the spatial correlation. The authors of [33] proposed the switching median filter (SM) algorithm, and a further optimized progressive switching filter (PSM) algorithm based on the SM algorithm is proposed in [34]. The minmax algorithm and the weighted median filter (WM) algorithm are proposed in [35,36], respectively. These algorithms have allowed useful explorations in improving the performance of median filters, but have their own limitations in practical applications because of the different interferences or noises that they deal with. The SM algorithm works well when dealing with low interference or noise conditions, but its performance gradually approaches that of the standard median filtering algorithm as the signal-to-noise ratio (SNR) of input images decreases. The PSM algorithm is a cyclic operation, which takes a long time to execute, and it has the requirement of estimating parameters in advance, which limits its real-time application. Although the minmax algorithm reduces the accumulation of propagation of error to a certain extent, the blurring of details is not well resolved. The WM algorithm reduces the loss of details by weighting, but its denoising performance also decreases.
Taking the EM algorithm as an example, the filtering method can be expressed as where x ij denotes each pixel of the image, i and j denote the row and column coordinates of the pixel, respectively, W[x ij ] denotes a window of points in the image x ij , med(W[x ij ]) means taking the median of all points in the window W[x ij ], y ij is the output of image x ij processed by the EM algorithm, N is the set of noise points and S is that of signal points. The EM algorithm shows advantages in both speed and performance. However, most of the improved median filtering algorithms, including the EM algorithm, are oriented to the identification of interference or noise points, rather than to the judgment of its ratio. In fact, due to the change in the ratio of co-channel interference, the existing median filtering algorithm cannot effectively remove them in high-ratio regions.

Model Construction
After using the improved Laplace operator to detect the co-channel interference and binarizing the original radar image I1000×3000 , the interference judgment matrix B1000×3000 is obtained. A 60-column row vector is used as a sliding window to traverse the entire radar image, and the tangential interference ratio is stored in turn to construct a tangential interference ratio model. The interference points are divided into high-ratio regions and low ones by setting a threshold. Among this, we complete the classification with the inter-class variance. The interference ratios obtained in previous steps are divided into two categories; the within-class variance of the two categories is calculated by the method of second-order cumulative moment. The largest inter-class variance between them is then obtained by class separability measurements as the division principle of high-ratio regions and low ones. For the interference points in low-ratio regions, the adaptive median filtering algorithm is used to remove the co-channel interference. For high-ratio regions, the adaptive median filtering algorithm is further improved. Starting from the defined minimum window size, we determine whether median value of the current sub-window is an interference point. If not, we take its gray value and assign it to the current pixel point; otherwise, we expand the window size to judge again. When the last ∆n pixels of the filter window are removed, this is the sub-window defined in this paper. In this process, ∆n is the difference between the actual tangential interference ratio in the region and the ratio threshold. A flowchart of the improved adaptive median filtering algorithm based on between-class variance can be represented by the following figure.
After obtaining the interference judgment matrix, the first step is to construct a tangential interference ratio model. By setting a threshold, the interference points are divided into high-ratio regions and low ones. The tangential interference ratio model is a graph drawn according to the statistical law of the interference ratio in the x-direction in the Cartesian coordinate system. A 1 × 50 sliding window is used to traverse each row of the original radar image matrix in turn; the ratio of interference points to the pixel blocks of the window is saved in turn and then drawn into a curve graph. In this paper, we take 4.4 n miles, 5.3 n miles and 6.1 n miles to the radar detection center as examples, and the tangential interference ratio model is expressed as follows.
The dotted line in magenta in Figure 2 represents the interference ratio everywhere along the tangential direction at a distance of 4.4 n miles. Similarly, the solid line in cyan and the dot dash line in blue represent the interference ratio at 5.3 n miles and 6.1 n miles, respectively.
In the second step, since the existing median filtering algorithm has a poor denoising effect on the dense regions of co-channel interference, this paper proposes a method based on between-class variance to find the tangential interference ratio threshold to distinguish high-ratio regions from low-ratio regions. Suppose that the tangential interference ratio obtained in the previous step is represented in ρ levels (0, 1, 2, . . ., ρ). Let ρ i denote the number of interference points at ratio level i, and D denote the total number of interference points [37]. The probability of occurrence of level i is given by Let the tangential ratio of interference points be divided into ρ L and ρ H according to the threshold ρ N . ρ L consists of interference points with a tangential interference ratio level of [0, . . ., ρ N ], and ρ H consists of interference points with a tangential interference ratio level of [ρ N + 1, . . ., ρ]. The innovation of this paper is mainly aimed at the processing of the points with a level of ρ H . Let P H (ρ N ) and P H (ρ N ) denote the cumulative probability. µ L (ρ N ) and µ H (ρ N ) represent the mean levels. σ 2 L (ρ N ) and σ 2 H (ρ N ) denote the variances of the classes ρ L and ρ H , respectively. The specific calculation rule is given as Figure 3.  The probability of occurrence of a interference point with a tangential interference ratio [0, . . ., ρ N ] can be given by Similarly, the probability of occurrence of a interference point with a tangential interference ratio [ρ N + 1, . . ., ρ] is We divide the tangential interference ratio into ρ L and ρ H , according to the knowledge of probability theory [38]: where µ L (ρ N ) and µ H (ρ N ) are the mean values of these two types of tangential interference ratio regions, respectively, and µ is the global mean value.
According to the concept of variance, the expression of between-class variance is According to Equation (13), we traverse each tangential interference ratio obtained in the first step in turn [39]. When σ 2 b (ρ N ) takes the maximum value max{σ 2 b (ρ N )}, the threshold of the high and low tangential interference ratio regions ρ L and ρ H is defined. The algorithm used to accomplish this task is as follows Algorithm 2.

Algorithm 2 Tangential interference ratio model construction
Step 1: Obtain first address of the original image and size of the image; Step 2: Open up a memory buffer to temporarily store the processing results and initialize it to 0; Step 3: Obtain the binarized image and scan it with a certain row vector window; save the ratio of interference points in the window; Step 4: Scan pixels in the binarized image through the loop statement; Step 5: Repeat step 4 until all pixels in the binarized original image are processed; using the maximum between-class variance method to find the threshold of the interference ratio, the obtained interference points are divided into high-ratio regions and low ones. By calculation, max{σ 2 b (ρ N )} = 9.8% is the tangential interference ratio threshold selected in this paper. The high and low tangential interference ratio regions segmented by this threshold are as follows.
The bright sparks in Figure 4a,b represent regions of high and low tangential interference ratio, respectively. In order to better understand Figure 4a,b, a part of the original image has been extracted and displayed on the left side of the corresponding picture. The left side of Figure 4a is the tangential interference points with interference ratio greater than 9.8%, while Figure 4b is the opposite.

Improvement of Median Filtering Algorithm
The method described in this section is optimized for regions of high tangential interference ratio in Section 3.3.
Due to the characteristics of high tangential interference ratio regions being greatly interfered with by the other interference in the window during the median filtering process, the adaptive median filtering in the window cannot fundamentally improve the interference removal effect. This section reduces the interference by adaptively changing the elements in the window. The core of the traditional median filtering algorithm is y ij = med(W[x ij ]) when a point is judged to be an interference or noise point [40]. This paper sorts the gray value elements in the selected window firstly to obtain Sort(W[x ij ]). Then, the last ∆n elements in the window are removed, where ∆n is the difference between the tangential interference ratio at this point and the ratio threshold ρ N . The new window obtained after this step is defined as W sub [x ij ] in this paper. This step draws most of the highbrightness interference in the window to the outside, thus greatly avoiding its influence on the filtering effect.
(a) (b) Figure 4. Comparison of high and low tangential interference ratio regions. (a) High tangential interference ratio regions, (b) low tangential interference ratio regions. Figure 5a is the processing result of the adaptive median filtering algorithm. The method in this paper can effectively suppress more co-channel interference, such as the part highlighted by the arrow in Figure 5a, thus showing that the improved adaptive median filtering algorithm has better performance.

Experiment and Simulation
In this paper, several commonly used image quality evaluation parameters are used to evaluate the advantages and disadvantages of the improved algorithm and several other filtering methods.
Let f (x, y) be the original radar image signal received by the radar receiver, andf (x, y) be the image signal processed by different algorithms. The Root Mean Square Error (RMSE) represents the average error between image pixels [41], which is defined as (14) where H and W represent the size of the image data.
The mesh surface graphs of the original radar image, the image processed by the adaptive median filtering algorithm and that processed by the algorithm in this paper are represented in Figure 6a-c. Compared with the adaptive median filtering method, the RMSE of the method in this paper is reduced from 3.8060 to 2.0561 by (14). The proposed algorithm removes some interference points in the filter window and optimizes the selection of the median, whereas the RMSE is sensitive to outliers in data. Therefore, through calculation, the proposed algorithm reduces the RMSE of the radar image by 1.7499, which improves the image quality.
Peak Signal to Noise Ratio (PSNR) is the most common and widely used image objective evaluation index; it is based on the error between corresponding pixels, i.e., an error-sensitive image quality evaluation [42]. The PSNR is expressed as follows: where PSNR is a measure of image quality. An increase in PSNR represents higher image quality. Compared with the adaptive median filtering method, the PSNR of the algorithm in this paper increased from 42.3 to 45 by (15). It can be seen that the proposed algorithm can better remove the co-channel interference in the radar image and improve the image quality.
When the Signal to Noise Ratio (SNR) is used to judge the quality of filtered images, the larger its value is, the more thorough the image filtering is. SNR is also a key indicator for evaluating image filtering quality, which can be given by By calculation, the radar images processed by the adaptive median filtering algorithm and the proposed algorithm have PSNRs of 10.8453 and 15.9134, respectively. The higher SNR indicates that the improvement of the proposed algorithm is more suitable for radar co-channel interference.
In order to prove the universality and effectiveness of the proposed algorithm in suppressing radar co-channel interference, the adaptive median filtering algorithm and the proposed algorithm are used to filter 20 radar images, and the RMSE, PSNR and SNR are calculated, as shown in Table 1. In order to better verify the effectiveness of method in this paper, we draw SNR line graphs of images processed by the adaptive median filtering algorithm and the proposed algorithm under different tangential interference ratios.
It can be seen from Table 1 and Figure 7 that the PSNR and SNR of the algorithm in this paper are both larger than those of the adaptive median filtering algorithm, which shows that the former is more thorough in filtering and has less impact on the ocean information of the original image. In addition, the effect of the proposed algorithm is more obvious as the tangential interference ratio increases.  Therefore, it can be concluded that the improvement of the adaptive median filtering algorithm in this paper is better than the algorithm before improvement, as it can deal with co-channel interference well and maintain the original ocean information. This improved image feature can be adapted to the later extraction of ocean information, especially the extraction of weak echo signals such as sea ice echoes.

Conclusions
The application of the median filtering algorithm and its various improved algorithms in the software removal of radar co-channel interference is studied. By improving the Laplace operator, a line detection template is formed, thereby enhancing the radial grayscale mutation of the radar image. In the filtering of interference points, the algorithm in this paper inherits the advantages of the adaptive median filtering algorithm. At the same time, the proposed algorithm establishes a tangential interference ratio model and determines the ratio threshold based on maximum between-class variance. For interference points in low-ratio regions, we replace the values of interference points with the median of current windows. Moreover, the median of the current window after removing the last ∆n elements is used to replace the value of interference points in high-ratio regions, so as to prevent the gray value of the dense interference point from interfering with the selection of the median value.
Results show that the proposed algorithm can effectively remove the co-channel interference in original radar image data. The RMSE of the image processed by the improved algorithm dropped from 3.8060 to 2.0561, while both PSNR and SNR increased. It is proven that the proposed algorithm has a more obvious effect on removing the co-channel interference, and has less influence on the useful signal of original images. The proposed algorithm can also be applied to the removal of interference or noise in other situations. For further research, due to the complexity of sea conditions, research based on the proposed algorithm and the echo characteristics of rain, snow and waves would be meaningful.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.