Validation and Evaluation of a Ship Echo-Based Array Phase Manifold Calibration Method for HF Surface Wave Radar DOA Estimation and Current Measurement

: Shore-based phased-array HF radars have been widely used for remotely sensing ocean surface current, wave, and wind around the world. However, phase uncertainties, especially phase distortions, in receiving elements signiﬁcantly degrade the performance of beam forming and direction-of-arrival (DOA) estimation for phased-array HF radar. To address this problem, the conventional array signal model is modiﬁed by adding a direction-based phase error matrix. Subsequently, an array phase manifold calibration method using antenna responses of incoming ship echoes is proposed. Later, an assessment on the proposed array calibration method is made based on the DOA estimations and current measurements that are obtained from the datasets that were collected with a multi-frequency HF (MHF) radar. MHF radar-estimated DOAs using three calibration strategies are compared with the ship directions that are provided by an Automatic Identiﬁcation System (AIS). Additionally, comparisons between the MHF radar-derived currents while using three calibration strategies and Acoustic Doppler Current Proﬁlers (ADCP)-measured currents are made. The results indicate that the proposed array calibration method is effective in DOA estimation and current measurement for phased-array HF radars, especially in the phase distortion situation.


Introduction
Up to now shore-based HF surface wave radars, which operate in the band from 3-30 MHz, have been extensively used for ocean surface current [1], wave [2], and wind measurement [3], as well as ship detection [4][5][6]. The capability of obtaining ocean parameters with a relatively high temporal resolution and a fairly extensive coverage makes HF radars suitable for oceanography research [7], hazard forecasting [8], coastal engineering [9], and maritime security [10,11]. Many types of HF radars are in operation worldwide [12][13][14]. Estimating the range and bearing from the backscattering patch to the receiving antennas is a necessary step for a HF radar to map ocean surface currents [15]. The distance from the backscattering patch to the radar is quite easy to estimate, whereas the estimation of target bearing is much more complicated. Direction-finding and beam-forming methods have been commonly adopted to estimate the direction-of-arrival (DOA) of the sea echo for HF radars [16].
For example, direction-finding techniques, such as multiple signal classification (MUSIC) algorithm, have been used by the Seasonde HF radar and the multi-frequency HF (MHF) radar to determine DOAs of sea surface current echoes. Kirincich et al. recently compared the performance of different direction-finding methods for HF radar current measurements [17]. In contrast, the WavE RAdar (WERA), which takes advantage of a long receiving array, adopts the beam-forming method to steer in the desired direction. However, in practical situations, both of the approaches are negatively affected by array uncertainties [18]. Without array calibration, the beam-forming method may result in an incorrect beam with a broadened beamwidth, and the direction-finding approach would yield a spurious direction with a poor angular resolution [19][20][21].
Unlike the use of array amplitude pattern in direction-finding HF radars [22], the performance of a phased-array HF radar significantly depends on the accuracy of the phase differences (or time delay) among receiving elements [23]. Cost-effective and real-time array calibration methods for phased-array HF radars have received significant attention in the past several years. For example, Fabrizio et al. developed an adaptive algorithm that estimates the receiver frequency response corrections for HF radar arrays [24]. Later, a maximum-likelihood (ML)-based phase error calibration method in which sea echoes are regarded as disparate sources was developed for the Ocean State Monitor and Analysis Radar (OSMAR) [25]. Afterwards, Fernandez et al. developed a phase error calibration method treating ship echoes as disparate sources [26]. Subsequently, Flores-Vidal et al. employed a number of unknown ship echoes to solve the phase errors from a cost function [27]. Recently, Chen et al. developed an approach that is based on the MUSIC algorithm (so-called MU method) for calibrating phase errors while using single-DOA sea echoes [28], and this method was compared with the ML-based method based on current measurements [18]. Almost all of the aforementioned methods only consider the phase uncertainty for each receiving element as a constant value that does not vary with direction. However, phase distortion in the receiving elements [29,30], which has a negative impact on the radar performance, is neglected in the aforementioned methods. This phase distortion problem was recently found in a small-aperture HF radar based on a circular array [31].
In order to address the phase distortion problem, the conventional array signal model is modified by adding a direction-based phase error matrix, and a new array phase manifold calibration method using measured phase responses of incoming ship echoes is proposed to compensate the array phase distortions. An assessment on the proposed array calibration method is made based on the DOA estimations and current measurements that are obtained from the datasets collected by a multi-frequency HF (MHF) radar, especially in the phase distortion case. The paper is organized, as follows: Section 2 contains the array signal model and a modified array signal model. The array phase manifold calibration method and DOA estimation method are also presented in this section. A description of the MHF radar and an example of array phase manifold measurement are given in Section 3. The assessment on the proposed calibration method based on DOA estimations and current measurements is provided in Section 4. The results are discussed in Section 5. Section 6 consists of a conclusion.

Conventional Calibration Model
Consider a narrowband signal s(t) with a DOA of θ observed by a receiving array consisting of M elements, in the absence of the mutual coupling, the signal at the output of the mth antenna element [28], z m , can be described by z m (t) = g m e jφ m e jωτ m s(t) + n m (t), where the subscript m denotes the mth antenna element, t represents time, g m = 1 + α m , α m and φ m are the gain and phase errors, respectively, ω denotes the radar operating radian frequency, τ m = (x m sin θ + y m cos θ) /c represents the time delay of radio wave, (x m , y m ) defines the position of the mth antenna element for a planar array, c represents the speed of light in free space, and n m (t) is the additive receiver noise. Let A m = e jωτ m , then the output vector Z(t) of M antenna elements can be expressed while using vector notation as where T is the transpose operator, Φ = diag [g 1 e jφ 1 , g 2 e jφ 2 , · · · , g M e jφ M ] represents the gain and phase error matrix, A = [A 1 , A 2 , · · · , A M ] T denotes the steering vector, and N (t) = [n 1 (t), n 2 (t), ..., n M (t)] T represents the additive noise vector.
In this model, the gain and phase errors are assumed to be constant for each antenna element. The model has been widely used for phased-array HF radars. Some aforementioned methods, e.g., [27,28], have been developed to estimate the gain and phase errors. For this model, it is straightforward to compensate the received signals if the gain and phase errors Φ have been estimated. Subsequently, the output vector of M antenna elements becomes From the output vector, the DOA of the radar returns can be estimated using a direction-finding approach [1]. Actually, the phase distortion exists in the receiving elements. Figure 1 shows an example of the phase values of A m measured from an MHF radar operating at 16.8 MHz. In Figure 1, the red lines, which are simulated using (3) from the receiving array configuration of MHF radar (see Figure 2), denote the ideal phase values against the signal direction. In contrast, the blue dots, which are measured from the known ship echoes that were collected with the MHF radar, represent the measured phase values. If there is no phase distortion, the correlation coefficient (CC) between the ideal phase values and the measured phase values should be 1. However, the CCs for antenna 3 and antenna 4 are 0.76 and 0.89, respectively. Obviously, there are some differences between the ideal and the measured phase values, and this difference is a function of direction.

Modified Calibration Model
The phase distortion leads the phase differences in Figure 1, and the gain and phase errors vary with direction in the presence of phase distortions. Consequently, the array signal model in (1) for a source signal coming from θ should be modified using direction-based gain and phase errors as z m (t) = g m (θ)e jφ m (θ) e jωτ m s(t) + n m (t).
Subsequently, the output vector of M antenna elements in (2) can be expressed as where the gain and phase error matrix Φ(θ) is a function of θ and Φ(θ)= diag [g 1 (θ)e jφ 1 (θ) , · · · , g M (θ)e jφ M (θ) ] . It can be found from (4) that, for a given θ, the value is constant. Let the first antenna element be a referenced element, then (4) could be expressed as and z 1 becomes where B 1 = 1. From (7) and (8), the output vector of M antenna elements in (5) becomes where B = [g 1 B 1 , g 2 B 2 , · · · , g M B M ] T denotes the array phase manifold or the new steering vector. The gain and phase errors in (5) are associated with direction, and the bearing of the incoming signal is unknown. Therefore, it is impossible to directly compensate the received signal to obtain the signals in (5). Alternatively, the term B(θ) can be treated as a steering vector and estimated while using the responses of receiving antenna array. If the phase of B(θ) is measured (like the blue dots in Figure 1), the B(θ) can be calculated and then used to estimate a more accurate DOA rather than using the term A(θ).

Array Phase Manifold Estimation
For a known ship echo with a relatively high signal-to-noise ratio (SNR), e.g., 20 dB, the estimate of array phase manifold, which is denoted as B m (θ), can be obtained from receiving signals as: The phase angle of B m (θ), which is denoted as ϕ m (θ), can be estimated as According to (10) and (11), a method for estimating the array phase manifold B(θ) while using antenna responses of incoming ship echoes (abbreviated as APM method) is proposed. The steps are as follows: 1. Use the MU method to compensate the phase errors of the antenna array. 2. Detect ships by applying a constant false alarm rate (CFAR)-based method [32] to the range-Doppler spectra. 3. Discard the detected ship echoes in the first-order clutter zone. The Doppler frequency of the first-order clutter zone is assumed to be in the band from −1.3 f B to −0.7 f B and from 0.7 f B to 1.3 f B ( f B denotes the Doppler Bragg frequency). 4. Obtain ship information, such as time, range, and radial velocity from both radar and AIS data. 5. Determine whether the radar echoes obtained in Step 3 are ship echoes using AIS data. 6. Calculate the phase data of the ship echoes by using Equations (10) and (11). From the AIS data the Doppler frequency of the individual ship, f ais , can be calculated. If the ship peak exists in the frequency range from f ais − 5 * ∆ f to f ais + 5 * ∆ f (∆ f denotes the Doppler frequency resolution), the individual ship is matched to the AIS data. 7. Apply a least-squares fitting to ship echoes to obtain the phase angle of B m (θ), then the array phase manifold B m (θ) is calculated .
Because the bearing distribution of incoming ship echoes is random, a least-squares fitting method based on a polynomial model is employed to estimate the complex value B m (θ) in Step 7. This polynomial model can be expressed as where x denotes the bearing of incoming signal, y represents B m (θ) in (11), and b is the order of the polynomial function (b = 30 in this work). For a small b (e.g., 4, 10, or 15), the fitted B m (θ), which is obtained using (12), can not match well with ship echoes. In contrast, for a big b (e.g., 50), it may need more ship echoes. It seems that b = 30 is appropriate to obtain B m (θ) from hundreds of ship echoes for the MHF radar. It can be seen from the steps that the approach takes few seconds while using a computer and needs accessional AIS antenna and receiver.

DOA Estimation for Ship Detection and Current Measurement
The MUSIC algorithm is used for HF radar DOA estimation and current measurement. For phase error calibration, the phase error of the received signals should be compensated first (i.e., Z in (3) is obtained). Subsequently, the covariance matrix R of the received signals can be obtained where H denotes the Hermitian transpose operation and E represents the mathematical expectation. The spatial spectrum, P, obtained from R in (13) while using the MUSIC algorithm can be expressed by (e.g., see [17]) where E N represents the noise subspace derived from the eigen-decomposition of the covariance matrix R. In contrast, the utility of the measured array phase manifold for the MUSIC algorithm is different. The term A in (14) should be replaced by the measured array phase manifold as: Eight eigenvalues satisfying from the eigen-decomposition of the covariance matrix. The biggest gradient of the eigenvalues (k 1 /k 2 , k 2 /k 3 , ..., k 7 /k 8 ) is calculated as the boundary between the signal subspace and noise subspace.

MHF Radar
The datasets that were selected in this work were collected with a MHF radar installed at Zhujiajian along the coast of the East China Sea from 28 August to 1 September 2009. Table 1 provides the details on the radar parameters. The transmitting antennas and receiving antennas are separated. The receiving array consisting of eight monopole antennas requires an occupying area of 54 × 12 m 2 , and the positions of the receiving antennas are shown in Figure 2. The steering vector of radar array can be theoretically calculated using the positions of the antennas. During the experiment, the MHF radar operated in a time division multiple sweep mode at 12.5 MHz, 16.8 MHz, and 19.1 MHz, with range resolutions of 5 km, 2 km, and 1 km, respectively. A 512-point fast Fourier transform (FFT) with a 75% overlap is applied to the raw radar data to obtain the Doppler spectra. Every covariance matrix consists of eight snaps. Therefore, 1408 chirps are involved in order to estimate the DOAs while using the MUSIC algorithm.

Array Manifold Measurement
There is a harbor near the Zhujiajian radar site, and lots of ships pass through the radar coverage everyday. Figure 3a,b show two range-Doppler spectra that were collected at 16.8 MHz and 19.1 MHz at 10:40 a.m. on 28 August 2009, respectively. Each range-Doppler spectrum contains seven peaks that are caused by ships. In an attempt to measure the array manifold using ship echoes, an AIS, which is an automatic tracking system that is commonly used on ships and by vessel traffic services (VTS) for identifying and locating vessels, was installed near the radar station during this experiment. The AIS system provided real-time information of ships (e.g., position and speed) within an approximate 100 km radius. In the meantime, the Zhujiajian MHF radar received the echoes of corresponding ships. Therefore, combining AIS and MHF radar data together, ship echoes, as well as the ship locations relative to the radar array, were collected. The ship echoes outside the first-order clutter zone, which may not be masked in sea clutter, are convenient to identify. A minimum SNR of 20 dB is required for array phase manifold measurement in this study. It can be found that array phase manifold is also a function of radio wavelength. Therefore, the array phase manifold needs to be measured at every radio frequency. Figure 4 shows an example of ships that meet the requirement to measure array phase manifold at 19.1 MHz. The distances between the selected ships and the radar are mostly less than 50 km.  Figure 5 shows the phase angles of the measured array phase manifold at 19.1 MHz. The red dashed lines, which are simulated using (6) from the receiving array configuration of the MHF radar as Figure 2, denote the ideal phase angles of B m (θ) against signal direction θ. The blue dots represent the measured phase angles ϕ m , which are obtained from ship echoes using Equations (10) and (11). The black lines are acquired by applying a least-squares fitting to blue dots. Consequently, the array phase manifold at operating frequency of 19.1 MHz is obtained using actual phase angles against signal directions relative to antenna boresight from −15 to 35 • (the black lines). Note that phase errors have been compensated in the figure. The fitting values (black solid lines) roughly coincide with theoretical values (red dashed lines). However, in some bearing range (as the green rectangle shown in Figure 4, approximately from 13 to 20 • ), the fitting values are not in agreement with the theoretical values. This disagreement is caused by the phase distortion of antenna array. In the phase distortion situation, the theoretical values would bring errors to DOA estimations. By the same way, array phase manifolds at operating frequencies of 12.5 MHz and 16.8 MHz are also measured using ship echoes, as shown in Figures 6 and 7. Figure 8a-c show the number of AIS points with ship DOAs. Figure 8d-f show the mean differences between the ideal phase and the fitting phase against direction at 12.5 MHz, 16.8 MHz, and 19.1 MHz, respectively. In the phase distortion range (as the black rectangle shown), the mean differences are larger than those in the other ranges.

Ship DOA Estimation
Ship DOAs are estimated from the MHF radar data using the MUSIC algorithm, which is an eigenstructure-based method, in order to evaluate the performance of the proposed calibration method. The ship bearings provided by the AIS are treated as 'true direction'. As shown in Figure 9, DOA estimation comparisons are made among three data calibration strategies: no calibration, phase error calibration, and array manifold calibration. No calibration means that the radar datasets have not been compensated before DOA estimation. Phase error calibration represents that phase errors of the radar have been compensated using a MU method before DOA estimation [28]. The MU method forms a cost function using the eigenstructure decomposition of single-DOA first-order sea echoes, and the phase errors in array are estimated by minimizing the cost function through an iterative procedure. For phase error calibration, the MU method and the ML approach have already been compared in the previous work [18], and that work indicates that the MU method is better than the ML approach for the MHF radar. Therefore, the MU method is adopted for comparison in this paper. Array manifold calibration means that the array phase manifold has been measured while using the APM method that was introduced in Section 2.3, and ship DOAs are estimated using the method presented in Section 2.4.   , 467, 1638, and 1855 ship echoes are, respectively, used for comparisons for the radar datasets at operating frequencies of 12.5 MHz, 16.8 MHz, and 19.1 MHz. The root mean square deviations (RMSDs) of ship DOA estimation using the no calibration approach at these three operating frequencies are all larger than 18 degrees, as shown in Table 2. In contrast, both calibration strategies using the MU method and the APM method significantly improve the DOA estimation performance compared with the no calibration approach. The results indicate that the APM method is the best among these three strategies with RMSDs around 1.5 degrees at these operating frequencies. For the 19.1 MHz datasets, significant phase distortion approximately occurs in the bearing range of radar array from 13 to 20 • , as shown in Figures 5 and 8f. The datasets at 12.5 MHz and 16.8 MHz have also been tested. In the radar bearing range from −7 to 10 • at 12.5 MHz and from 3 to 15 • at 16.8 MHz, the phase distortion is significant. As shown in Table 3 and Figure 10, in the radar bearing range where the phase distortion is significant, the ship DOA measurement accuracy using the phase error calibration method drops. The APM method shows better performance than the phase error calibration method for the phase distortion cases. The RMSDs are reduced from 2.7 degrees to 1.8 degrees at 12.5 MHz, from 2.2 degrees to 1.6 degrees at 16.8 MHz, and from 1.9 degrees to 1.3 degrees at 19.1 MHz. In the case that phase distortion is significant, the RMSDs have been reduced around 30% while using the APM method compared with the phase error calibration method. The slopes and biases of ship DOA measurements at 12.5 MHz, 16.8 MHz, and 19.1 MHz are also shown in Figure 10. When comparing with the phase error calibration method, the slopes of the ship DOA measurements using the APM method are closer to 1, and the biases are smaller. This also indicates that the APM method is better than the phase error calibration method.

Current Measurements
The MHF radar was employed to provide current maps over a portion of the East China Sea near Zhejiang Province. Currents in this area are complex, and the radial current v R obtained from the MHF radar is a function of radar look direction θ R , i.e. v R = vcos(θ C − θ R ), where v is the current velocity and θ C denotes the current direction. Consequently, making a comparison of radar-measured radial current against 'ground-truth' data is a convincing way to evaluate calibration methods that significantly affect the estimation of radar look direction θ R . The radial currents obtained from the MHF radar data using three calibration strategies (no calibration, phase error calibration, and array manifold calibration) are compared with Acoustic Doppler Current Profilers (ADCP) current measurements that are considered as 'ground-truth' data.
A bottom-mounted upward-looking ADCP was deployed at the cell that is 10 km from the radar and −10 • from the radar boresight. The ADCP provided current vector profiles at depths from approximately 3 m to 20 m at a resolution of 1 m and an interval of 20 min. The ADCP-measured current vectors at the depth of 3 m are converted to radial currents for comparisons.
Radar-derived radial currents are estimated from all eight receiving elements data while using the MUSIC algorithm every 10 minutes. Details on the MHF radar current measurements can be found in [1]. A first-order signal SNR threshold of 10 dB and a MUSIC spectrum SNR threshold of 6 dB are configured for quality control. Figure 11a-c show the radial current comparisons at operating frequencies of 12.5 MHz, 16.8 MHz, and 19.1 MHz, respectively. CCs and RMSDs are reported in Table 4. When compared with the no calibration approach, both phase error calibration and APM methods significantly improve the current measurement performance. For example, at 16.8 MHz, the current measurements without array calibration are poor with a CC of 0.40 and a RMSD of 22.4 cm/s, whereas current measurements obtained by using the phase error calibration and the APM methods are consistent with ADCP-measured currents, with CCs of 0.82 and 0.86, respectively. The CCs and RMSDs indicate that the measurements obtained using the APM method are better than those using the phase error calibration method. Additionally, the comparisons at operating frequencies of 12.5 MHz and 19.1 MHz indicate that the APM method is the best.

Current Measurements Using Four Receiving Elements
For a special antenna configuration, radial currents are obtained from four-element MHF radar datasets. The four receiving elements are antenna 2, antenna 3, antenna 4, and antenna 5, which are deployed as the configuration that is shown in Figure 2. The experiment is conducted in this way, because the four-element configuration has been employed in some phase-mode direction-finding radars, e.g., the compact WERA system that was developed by the Heltz company and the High Frequency Doppler Radio scatterometers (HFDR) in direction-finding mode developed by University of Hawaii. Radial currents which are obtained from the MHF radar datasets at the location of ADCP using three calibration strategies are compared with ADCP measurements. Figure 12a-

Discussion
HF radar current measurement and target detection performance is associated with DOA estimation capability, which strongly depends on the accuracy of array manifold. Both the proposed calibration technique (APM method) and traditional phase error calibration method (MU method) show significant improvements on HF radar DOA estimation and current measurement comparing with the no calibration strategy. In the case that phase distortion is significant, the RMSDs have been reduced while using the APM method as compared with the MU method.
For current comparison, the ADCP was deployed at the cell, which is −10 • from the radar boresight. At this bearing the phase distortions of the radar at these three operating frequencies are not significant (see Table 3). At the bearing of −10 • , the steering vector for the MU method A(−10 • ) and the steering vector for the APM method B(−10 • ) are approximately the same. As a consequence, the DOA performances at the bearing of −10 • using the APM method and the phase error method are nearly the same. This is why the APM method slightly, but not significantly, improves the current measurement performance compared with the MU method for these datasets.
Another interesting work is that the APM method has been applied to four-element HF radar datasets. As reported in Table 5, the current measurements that were obtained from the four-element HF radar datasets using the APM method are as good as those obtained from eight-element MHF radar datasets using the MU method. This is extremely useful for compact HF radars, since they are more convenient than large array radar to deploy.

Conclusions
The receiving array signal model for a phased-array HF radar has been investigated and modified. Phase distortions have been observed in the MHF radar receiving elements. The APM method based on measuring array manifold HF radar ship echoes as well as AIS data has been proposed. The performance of the APM method has been compared with the no calibration strategy and the phase error calibration method.The DOA estimations and current measurements using the MUSIC algorithm show a reasonable improvement after taking the phase distortion into consideration and using the APM method.
Array calibration is critical for HF surface wave radar applications, and it can significantly improve azimuth estimation performance as well as current measurements. The phase error calibration, which has been studied in the past few decades, is not sufficient for the MHF radar in the phase distortion situation. The APM method provides a feasible way to measure the array manifold and calibrate the array distortion. This method will be validated in other radar sites in the future.