Improved Amplitude-Phase Calibration Method of Nonlinear Array for Wide-Beam High-Frequency Surface Wave Radar

Abstract

The amplitude and phase errors in array elements lead to uncertainties in the estimation of arrival angles (DOA), which, in turn, affect the accuracy of current measurements for high-frequency surface wave radar (HFSWR). To address this issue, this paper proposes a passive amplitude-phase self-calibration method for a wide-beam HFSWR system with a nonlinear array. This method utilizes the aggregation of the amplitude ratio among array elements to screen reception matrices of single DOA sources. Based on the differences in reception matrices, the amplitude error is calibrated. Moreover, the cost function is calculated using the multiple signal classification (MUSIC) algorithm, and the initial phase error is first obtained after triangular array dimensionality reduction. Then, the phase error is further calibrated through quadratic-form iterative optimization. This method has been validated through simulations and real measurements. An approximately 4-day dataset obtained via HFSWR is reanalyzed in this paper. After data calibration, the radar-estimated currents were in good agreement with the buoy-measured results, with a root mean square error (RMSE) as low as 0.06 m/s and a correlation coefficient (CC) of up to 0.88. The results indicate that this method is suitable for the amplitude and phase error calibration of nonlinear arrays in wide-beam HFSWR systems.

1. Introduction

Since the pioneering work of electromagnetic scattering from the ocean surface, high-frequency surface wave radars (HFSWR) located on the coast have been widely applied in marine environment monitoring over the world for wave, wind, and especially current measurements [1,2,3]. The resonance scattering mechanism between electromagnetic waves and the sea surface results in two discrete peaks in the Doppler spectrum [4,5]. The radial current maps can be obtained from these backscattering echoes [6].

The accuracies of the obtained current are influenced by many factors such as radio wavelength, power spectrum frequency resolution, and the direction of arrival (DOA) estimation performance of the first-order ocean echoes [7,8]. Radio wavelength and power spectrum frequency resolution can be easily determined via radar waveform parameters, while the DOA is always estimated by the multiple signal classification (MUSIC) algorithm and beamforming algorithm [9,10]. The algorithms have high resolution performance under the premise that the echo signal satisfies the hypothesis of ideal model. However, the array performances in practical applications will be significantly reduced for the array errors resulting from position, mutual coupling effect, or channel amplitude and phase errors [11]. Without effective calibration, a HF radar system not only estimates spurious directions but also provides a poor azimuth resolution for wave, wind, and especially current measurements. As a consequence, the array calibration algorithm has become increasingly important for HF radar ocean surface remote sensing.

A common method for array calibration can effectively correct the array manifold by placing a responder on the shore or ship as an auxiliary calibration source [12,13]. However, this method struggles to provide continuous real-time calibration values, and the installation and maintenance of the transponder are also relatively expensive in terms of economic costs. In recent years, the passive amplitude and phase calibration has been widely applied. Under the passive mode, without the installation of auxiliary equipment or the acquisition of prior information, self-calibration of amplitude and phase errors can be achieved by directly processing radar echo information. Meanwhile, the calibration will be approximately real-time while the ocean echo is stable. Friedlander and Weiss proposed a classical method for amplitude-phase error self-calibration by utilizing the orthogonality of noise subspace and signal subspace [14,15]. Subsequently, Solomon extended this method and took ionized meteor trajectories as single sources for passive amplitude and phase calibration [16]. In 2006, Wu Xiongbin utilized the shift-invariant antenna pairs to select the single sources in ocean echoes, and calibrated it through the maximum likelihood (ML) cost function [17]. The search for single DOA signals using shift-invariant antenna pairs is a relatively reliable method, while the maximum likelihood cost function requires multi-dimensional search for both direction of arrival and amplitude-phase error, which leads to the problem of overestimation of the parameters. In addition, not all arrays of HF radar have shift-invariant antenna pairs [18,19]. Similar to the ML method, Chen et al. proposed another appoach based on the multiple signal classification (MU) algotithm for amplitude and phase calibration using the reception matrices of single DOA sources [20]. Later, Zhao Chen validated and compared the ML method and the MU method, and the results showed that both methods can effectively improve the accuracy of ocean current estimation [21].

In summary, the passive amplitude and phase calibration method can provide better application prospects for estimating ocean current parameters for HF radar. Based on the perspective of orthogonality between noise subspace and signal subspace, this paper proposes a method for amplitude and phase calibration of nonlinear arrays in wide-beam HF radar. The reception matrices of single DOA sources are screened by the aggregation of the amplitude ratio between array elements. Combining triangular array dimensionality reduction with the MUSIC method, stable and accurate estimation of the initial phase values was achieved. To validate this method, the datasets collected with a OSMAR071G radar deployed at Pearl River Estuary in April 2021 were reanalyzed.

The article is organized as follows. Section 2 provides a description of the array signal model. The amplitude and phase calibration method is also given in this section. In Section 3, the method is tested using a simulated dataset. To validate the proposed method, four-day comparisons between radar-retrieved and buoy-measured current velocities are made in Section 4. Conclusions and perspectives on future work are presented in Section 5.

2. Array Signal Model and Amplitude and Phase Calibration Method

2.1. Array Signal Model

Assuming the array HF radar has M array elements, and the coordinates of array element i are ($x_i$, $y_i$) ($i=1,2\dots M$), $x_1=0$, $y_1=0$. After two fast Fourier transforms (FFTs), the echo signal model on channel i can be simply described as:

$$z_i(m,o)=g_i{e}^{\left(j{\phi }_i\right)}{\sum }_{k=1}^K{S}_k(m,o)a_i\left({\theta }_k\right)+n_i(m,o)$$

(1)

where m is the discrete frequency on the distance dimension; o is the Doppler frequency shift; $g_i$ and $e^{j{\phi }_i}$ are the amplitude and phase errors, respectively; $S_k$ represents a two-dimensional echo signal of an arrival angle of ${\theta }_k$; and $n_i$ is the additive noise. The steering vector $a\left({\theta }_k\right)$ can be represented by

$$a\left({\theta }_k\right)={[1a_2\left({\theta }_k\right)\cdots a_M\left({\theta }_k\right)]}^T,$$

(2)

where $a_i\left({\theta }_k\right)=e^{j\frac{2\pi }{\lambda }}(x_{M}sin{\theta }_k+y_{M}cos{\theta }_k)$, $\lambda$ is the radio wavelength.

The response of the receiving array is

$$X=G\Phi AS+N,$$

(3)

where $X=[Z_1{Z}_2\cdots Z_M){]}^T$, which represents an array snapshot during the coherent accumulation period; $G=diag\left({[1g_1\cdots g_M]}^T\right)$, $\Phi =diag{\left([1e^{j{\phi }_2}\cdots e^{j{\phi }_M})\right]}^T)$, which are amplitude and phase error matrix, respectively; $S={[S_1{S}_2\cdots S_k]}^T$, is the signal vector; $N=[n_1{n}_2\cdots n_M){]}^T$, is the array additive noise; and A is the theoretical manifold, $A=[a\left({\theta }_1\right)a\left({\theta }_2\right)\cdots a\left({\theta }_k\right)]$

2.2. Amplitude and Phase Calibration Method

According to Equation (1), when K = 1, the two-dimensional echo is a single direction of arrival signal:

$$z_i(m,o)=g_i{e}^{j{\phi }_i}S(m,o)a_i\left(\theta \right)+n_i(m,o).$$

(4)

If noise is ignored (i.e., $n_i(m,o)=0$), the above equation is simplified as

$$z_i(m,o)=g_i{e}^{j{\phi }_i}S(m,o)a_i\left(\theta \right).$$

(5)

Due to $g_1=e^{j{\phi }_1}=a_1\left(\theta \right)=1$,

$$z_1(m,o)=g_1{e}^{j{\phi }_1}S(m,o)a_1\left(\theta \right)=S(m,o).$$

(6)

Therefore, for a specific single direction of arrival signal, the amplitude error of array element i is

$$g_i=\frac{\left|z_i\right|}{\left|z_1\right|}=\left|\frac{g_i{e}^{j{\phi }_i}s(m,o)a_i\left(\theta \right)}{s(m,o)}\right|=\left|g_i{e}^{j{\phi }_i}{a}_i\left(\theta \right)\right|.$$

(7)

In fact, noise always exists, which results in different amplitude errors for different single direction of arrival signals. Taking the statistical average of the amplitude errors of multiple single arrival angles can further improve accuracy:

$$\widetilde{g_i}=\frac{\sqrt{{\sum }_l{\sum }_t{\left|z_i(l,t)\right|}^2}}{\sqrt{{\sum }_l{\sum }_t{\left|z_1(l,t).\right|}^2}}$$

(8)

After amplitude error correction, the array covariance matrix corresponding to any single arrival signal can be expressed as

$$R=E\left[X\left(t\right)X{\left(t\right)}^H\right]=\Phi a\left({\theta }_l\right){\sigma }_l^2{a}^H\left({\theta }_l\right){\Phi }^H+{\sigma }_n^{2}I,$$

(9)

where ${(.)}^H$ is the conjugate transposition; ${\sigma }_l^2$ and ${\sigma }_n^2$ are the signal power and noise power, respectively; and $X\left(t\right)$ is the recieved matrices corresponding to signal l. Decompose on R:

$$R={\sum }_{m=1}^M{\lambda }_m{\mu }_m{\mu }_m^H=E_s\left(l\right){\Lambda }_s{E}_s^H\left(l\right)+E_n\left(l\right){\Lambda }_n{E}_n^H\left(l\right),$$

(10)

where $E_n^H\left(l\right)$ is the noise subspace.

Due to $E_n^H\left(l\right)[\Phi a\left({\theta }_l\right)]=0$, the MUSIC sum cost function based on the orthogonal principle of noise subspace under single direction of arrival conditions is [22]

$$P=arg{\displaystyle \underset{\Phi ,{\theta }_l}{min}}{{\sum }_{l=1}^L∥E_n^H\left(l\right)\Phi a\left({\theta }_l\right)∥}^2.$$

(11)

Based on the eigenvalue structure of the cost function P, sufficient reception matrices of single DOA sources can be selected to ensure the universality and anisotropy of the arrival angle distribution in the case of uncertain single DOA sources.

When $K\ne 1$,

$$\frac{\left|z_i\right|}{\left|z_1\right|}=\left|\frac{g_i{e}^{j{\phi }_i}{\sum }_{k=1}^K{S}_k(m,o)a_i\left({\theta }_k\right)}{{\sum }_{k=1}^K{S}_k(m,o)}\right|.$$

(12)

Compare (12) with (7), for single DOA sources, the amplitude ratio between two antennae is exactly the amplitude error $g_i$ if noise is ignored, while it will concentrate around the amplitude error if noise is considered. However, for multi-DOA sources, the ratio does not show aggregation an effect since the modulus of ratio is random for the uncertainty of the arrival angle.

The aggregation effect mentioned above can be used for the selection of single arrival signals. Ideally, only two antennae are needed to find the single direction of arrival signal. However, the signal-to-noise ratio (SNR) for HF radar is usually around 20 $\mathrm{dB}$ or even lower. Therefore, if only two antennae are used, the above clustering effect will be weakened. Here, multiple antennae for joint estimation are considered.

A method is proposed by which to construct and screen reception matrices of single DOA sources in this article, utilizing the aggregation effect for multiple antennae mentioned above. Without considering phase information, the single DOA sources can be obtained solely through amplitude information, which means that there are no specific requirements on the antenna array layout. Figure 1 shows the algorithm flow. The steps are as follows.

Step 1: Divide sampling frames into multiple segments.

Step 2: Build a reception matrix via short-rime Fourier transform (STFT) for each segment. The frequency of the reception matrix corresponds to a two-dimensional frequency point after two FFTs. Set the number of receiving matrices as H.

Step 3: Construct a multidimensional spatial point set. For the reception matrix, divide each row of the matrix by the reference element row (generally the first row) and then take the modulus to form a new matrix. Average the matrix by rows to form a column vector of $M*1$. Combine all vectors into a M*H matrix. Select C (C > 2) rows with high-antenna SNR; then, a set of C-dimensional spatial points of size H will be formed for all the reception matrices.

Step 4: Find the single DOA sources by local clustering. Calculate the sum of squares of the distances from $p_i$ to the points which are the nearest $0.2*H$ points, denoted as $sum\_p_i$. When taking the minimum value of $sum\_p_i$, the $p_i$ is the local density center. $p_i$ and its surrounding $0.2*H$ points are identified as a single arrival.

To approximate the true phase error and the corresponding single arrival signal azimuth, it is necessary to optimize the minimum value of Equation (11). Based on the maximum likelihood cost function of the snapshot, the shift-invariant antenna pairs can be used to accurately determine whether a certain snapshot of the two-dimensional echo is a single arrival. However, due to the eigenvalue structure of the cost function P, the single arrival signal is not a specific snapshot but a reception matrix composed of continuous snapshots; and for L reception matrices, if the arrival angle distribution is limited to a too small interval, it will lead to ambiguity in the calibration values. Therefore, in the case of uncertain single DOA sources, sufficient reception matrices should be screened to ensure the universality and anisotropy of the arrival angle distribution.

In addition, to find a more accurate initial value of phase error, a triangular array dimensionality reduction and initial value-finding method based on the sum MUSIC cost function was designed. Figure 2 shows the algorithm flow. By constructing triangular subarrays and gradually adding elements, the complexity of the initial value estimation in Equation (11) will be reduced. Moreover, whether a signal is strictly a single arrival or not is not a necessary condition based on the idea of least squares. The setailed steps are as follows.

Step 1: Construct a nonlinear subarray (triangular subarray) with phase errors of $({\phi }_1,{\phi }_2,{\phi }_3)$. ${\phi }_1$ is the reference element phase error (typically set to 0). ${\phi }_2$ and ${\phi }_3$ are ergodically searched with a fixed step through $[-\pi ,\pi ]$ [23].

Step 2: Search for the arrival angle $\theta$ as well as the corresponding MUSIC spectral value.

Step 3: Obtain the initial phase error values of the array $({\phi }_1,{\phi }_2,{\phi }_3)$ and the arrival angles of each single source via Equation (11).

Step 4: Add the phase error of an adjacent array element with a fixed step from $[-\pi ,\pi ]$ to the array $({\phi }_1,{\phi }_2,{\phi }_3)$ and substitute the phase error values of the newly generated subarray determined via each traversal into Equation (11). Determine the phase error of the new array element by minimizing the cost function. Add array elements one by one until all the initial errors are determined.

Step 5: Calculate the optimal solution of phase error after quadratic-form iteration via Equation (11).

3. Simulation

Figure 3 shows a schematic diagram for an antenna, which has an eight-element dual nonlinear array. The horizontal spacing $d_1$ of array elements for the first row is set to be 18 m, for the second row, $d_2$, it is 36 m, and the vertical spacing $d_3$ between the two rows is 15 m.

Assuming the operating frequency of the radar is 8 MHz. Simulate the reception matrix corresponding to the two-dimensional echo spectral points using an 8 * 32 complex matrix. Amplitude errors from $g_1$ to $g_8$ are [0, 1.2, 2.85, −3, 0.92, 3.1, −2,5], respectively, and phase errors ranging from ${\phi }_1$ to ${\phi }_8$ are [$0^{°}$, ${25}^{°}$, ${37.2}^{°}$, $-{78}^{°}$, ${120.5}^{°}$, $31{}^{°}$, $-{50}^{°}$, ${77}^{°}$]. The flow pattern resolution of the spatial array used for estimating azimuth is $1^{°}$, the number of two-dimensional ocean echo signals is 2000, the arrival angle is randomly distributed between 1 and 4, and the single source ratio is one third. The triangular subarray used for estimating the initial value of phase error is 1-4-8, and the subsequent expanded antenna sequence is 2-3-5-6-7. Perform 50 Monte Carlo experiments using the method described in 2.2 under the SNR of 0 dB ∼ 50 dB and take the average of the calculated results as the amplitude and phase error result.

The accuracies of estimating single DOA sources when selecting different numbers of array elements under different signal-to-noise ratio conditions are shown in Figure 4. To maintain high accuracy, a sufficiently high signal-to-noise ratio is required for selecting single arrival sources from fewer array elements, while at lower signal-to-noise ratios, a large number of the array elements is needed. However, an increase in the number of antennae will lead to an increase in the computational complexity of the algorithm, and when the number of antennae reaches a certain threshold beyond which the improvement in accuracy is no longer significant. Therefore, practically, it is common to select a small number of array elements (≥3) with high SNR to estimate single DOA sources.

Figure 5 and Figure 6 show the trends of the estimation deviations of the amplitude error and phase error under different SNR conditions. As the signal-to-noise ratio increases, the deviation of amplitude error estimation decreases sharply. While the deviation of phase error estimation is not significantly reduced, it still shows a downward trend. Moreover, when the SNR is higher than 10 dB, the phase error estimation deviation is <1° and the amplitude error is ≤0.7 dB.

To further evaluate the performance of the proposed method, the impact of antenna position errors were taken into account. In Figure 3, the coordinates of antenna 1 changed from (−27, −15) to (−27.2, −15), the coordinates of antenna 4 changed from (−36, 0) to (−36.2, 0), and the coordinates of antenna 8 changed from (−63, −15) to (−63.1, −15). Although these errors are often neglected in engineering when the deviations are small, they still have a certain impact on amplitude and phase calibration.

Figure 7 and Figure 8 show the trends of estimation deviations for amplitude and phase errors under different signal-to-noise ratio (SNR) conditions for both the existing maximum likelihood (ML) method and the method proposed in this paper. Clearly, both the ML method and the proposed method demonstrate good calibration performance when the SNR is high (>20 dB). However, when the SNR is low, the ML method exhibits greater fluctuations, resulting in unstable calibration performance. This is because the ML method requires specific array layouts and the determination of the (or more) single arrival signals when searching for the initial values, which imposes certain limitations.

4. Validation

4.1. Equipment and Setup

From 30 April to 3 May 2021, an experiment was carried out to validate the amplitude-phase calibration method at Pearl River Estuary via a OSMAR071G radar. To realize monostatic operation, the radar employs frequency-modulated interrupted continuous wave (FMICW) chirps for radio energy radiation and uses a nonlinear array-receiving antenna, as shown in Figure 3, to receive ocean echoes. The radar waveform parameters are listed in

Table 1. Radar waveform parameters.

Parameters	Values
Operating frequency	8 MHz
Sweep bandwidth	30 kHz
Sweep period	0.6528 s
Arrary normal (based on due north clockwise)	170°
Range resolution	5 km
Angular resolution	1.5°
Detection range	250 km
Time resolution for current measurement	10 min

.

By measuring the projection of surface currents at any point on the sea surface towards the radar direction, radial currents can be obtained. A wave buoy was deployed about 60 km from the OSMAR071G. The positions of the radar and buoy are shown in Figure 9. The buoy-derived current measurements are presented as “sea-truth”. For convenience, both the buoy and radar produce current velocity every 10 min for comparison.

4.2. Observations and Discussion

We performed amplitude and phase calibration on echo data based on the method proposed in this article. Figure 10 and Figure 11 show the amplitude and phase calibration factors, respectively. Due to the weak signal-to-noise ratio (SNR) of the echo information from the single arrival angle used for calibration (below 12 dB), the datasets from 19:00 on 1 May to 22:00 on 1 May are not considered here.

A total of 200 single sources were selected at 14:30 on 1 May for amplitude and phase calibration via the method proposed in this article; then, we performed a MUSIC estimation on these 200 calibrated sources. Figure 8b and Figure 12a show the uncalibrated and calibrated MUSIC estimation spectra, respectively. There are many pseudo peaks in the uncalibrated MUSIC spectrum. It will be difficult to obtain the correct DOA since the spectrum has low SNR and wide beams, while the calibrated MUSIC spectrum, with high SNR, generally only has one peak and needle-like narrow beams, which makes it easy to determine the DOA.

Figure 13 shows the histogram statistics of the MUSIC-estimated arrival angles for the 200 single sources mentioned above after amplitude and phase calibration. Obviously, the arrival angle distribution of a single source is within the range of [−60, 60] centered on the array normal, which is also the radiation field angle range of the electromagnetic waves for the radar deployed at Wanshan. This demonstrates the rationality of the arrival angle distribution of a single source and further verifies the accuracy of amplitude and phase error estimation.

Comparisons of current velocity between the RADAR and the wave buoy have been undertaken. The comparison of radial flow fields of ocean currents with and without calibration at 9:10 am on 3 May is displayed in Figure 14. Due to the post-processing involved in ocean current inversion, the flow fields appear relatively smooth and well-defined under both calibrated and uncalibrated conditions. The radial flow field obtained from calibrated data gradually moves away from the radar as the longitude decreases, while the radial flow field obtained from uncalibrated data exhibits a very chaotic distribution in terms of flow velocity and direction. The time series are shown in Figure 15, including the buoy-measured results (blue-plus sign), the Radar-estimated results using the uncalibrated data (yellow-circle), and the Radar-estimated results using the data calibrated via the method proposed in this article (red-asterisk). The statistical characteristics are summarized in

Table 2. Statistical properties resulting from the comparisons of current velocities with the wave buoy.

	Uncalibrated	Calibrated
MD (m/s)	0.17	0.05
RMSE (m/s)	0.21	0.06
CC	−0.36	0.88

. Compared with the uncalibrated results, the current velocities extracted from the RADAR using the calibrated data agree better with the buoy-derived current measurements, with a CC of 0.88, an RMSE of 0.06 m/s and a mean difference (MD) of 0.05 m/s. The scatter plots of the uncalibrated and calibrated results are shown in Figure 16a and Figure 16b, respectively.

The signal-to-noise ratio (SNR) time series for the radar is shown in Figure 17a. The SNR varied between 10 and 40 dB. The relationship between the RMSE of the current measurements and SNR is also shown in Figure 17b. Obviously, the RMSE of the current velocities is less than 0.14 m/s and decreases with the increase in SNR.

5. Conclusions

An improved amplitude-phase self-calibration nonlinear array method based on a wide-beam HF radar system is proposed in this paper. As well as a simulation validation, a nearly four-day observation and intercomparison based on the OSMAR071G radar and buoy installed at Pearl River Estuary along the coast of Shen Zhen province in China has been analyzed in detail. The results indicate that the radar-derived measurements, after calibrating the echo data using the proposed method, have good agreement with the buoy-measured current velocity. To further validate the effectiveness of this method, we will consider the datasets under different sea conditions.