Mapping wind by the first-order Bragg scattering of broad-beam high-frequency radar

1. Introduction

High-frequency (HF) radar, as a shore-based remote sensing system, can monitor ocean state with the advantages of large-area, real-time and all-weather operations. The current and wave mapping measurement by HF radar are now a well-accepted technology (Zhou et al., 2015). However, mapping wind with HF radar is still a challenge. Ocean surface wind is a rich source of renewable energy due to its intensity and steadiness. Wind field extraction over the ocean is helpful for wind energy detection. In addition, offshore wind field information can help a wide variety of coastal and marine activities, such as sailing, fishing and surfing. Therefore, obtaining sea surface wind information from HF radar will make a meaningful contribution to wind energy detection, marine activities and oceanographic observation networks.

HF radar’s wind direction is estimated from the ratio of first-order Bragg scattering assuming that the Bragg wave is wind-driven and aligned with the wind direction. The reliability of this assumption depends largely on the selection of radar frequency, which corresponds to the Bragg wave frequency. For high radar frequency (short Bragg wave), the Bragg waves respond rapidly to local wind excitation and also decay rapidly when the winds cease. The short Bragg waves have the disadvantage that they will saturate rapidly and thus offer no information on the magnitude of the excitation (Dexter and Theodoridis, 1982). For low radar frequency (long Bragg wave), the Bragg waves saturate at a higher wind speed, but they respond slowly to local wind excitation and are very likely to be contaminated by swell. Because first-order scattering has a higher signal-to-noise ratio (SNR) than second-order scattering, the first-order based wind direction estimates have larger distance scales than the second-order based estimation. At present, the principle of wind direction estimation has been widely accepted, and its uncertainty lies in the determination of the wind spreading function (WSF) (Huang et al., 2004), and the left-right ambiguity relative to the radar beam. In general, generation of an HF radar wind direction map is relatively simple and feasible. Many researchers have obtained satisfactory wind direction estimation results by different HF radar systems (Heron and Rose, 1986; Fernandez et al., 1997; Wyatt et al., 2006; Wyatt, 2018; Huang et al., 2004).

Although inversion to obtain HF radar wind speed is more difficult than wind direction, some potential parameters of the HF radar Doppler spectrum that could be used have been proposed. The proposed parameters include (1) the ratio of the second-order power to the Bragg line power or the significant wave height (Barrick et al., 1974; Zeng et al., 2016), (2) the width of the spectral power surrounding the strongest Bragg line at a point 10 dB down from the peak (Stewart and Barnum, 1975), (3) the frequency position of the second-order peak (Green et al., 2009), and (4) the first-order Bragg peak (Shen et al., 2012; Kirincich, 2016). These methods can be mainly categorized into two types. The first type is that estimating the wind speed by the second-order scattering based on the second-order radar cross-section equation (second-order theory) (Barrick, 1972b) and the wind-wave relationship. The method has the advantage of not requiring additional calibration and having a wide unsaturated range. However, second-order scattering has a low SNR and is susceptible to external noise and spatial aliasing, which results in a small detection area. The second type involves estimation of wind speed from the first-order scattering only. This type includes using the relationship between the first-order spectrum power (FSP) and the energy in the directional wave spectrum at the Bragg wavelength and radar look direction (Barrick, 1972a) or using the first-order Bragg ratio (Shen et al., 2012). By using first-order scattering, these methods have some of the same features as wind direction estimation such as large detection area and rapid saturation relative to the second-order method. Because there is no calibration of the radar signal, FSP based wind speed extraction requires additional calibration from other wind detecting instruments. Maresca and Barnum (1982) concluded that the 10 dB width is not a good estimator of wind speed, and the wind speed estimates are better obtained from the second-order method by a prior wind-wave model. However, there are other restrictions of the second-order method in practical application besides the low SNR and short detection distance. These restrictions mainly come from the HF radar system technology and the details are described below.

HF radar industry offers two main genres for the technology: the phased-array genre and the crossed-loop genre. The phased-array genre has the benefits of narrow-beam and good temporal and spatial resolution, and the drawbacks of large size and high installation or maintenance costs (Heron, 2015). The crossed-loop genre has the merits of small antenna footprint and good economy, and the defects of broad-beam. The small antenna footprint facilitates a larger number of installations of crossed-loop systems around the world. Both genres can map wind direction well. However, for wind speed mapping, the two genres are very different. The phased-array genre can use both the second-order and first-order method to map wind speed. In calculating wind speed by the phased-array HF radar, the second-order method is widely used (Huang et al., 2002). For the crossed-loop HF radar, its second-order scattering cannot be extracted in each cell and it only provides homogeneous wind speed estimation in azimuth at a given range by the second-order scattering method. Hence to determine the wind speed field using this type of HF radar requires a first-order method.

The frequency shift of first-order scattering is fully used to extract surface current velocity. However, the intensity of the first-order scattering is typically neglected because of its saturation characteristic. In 2012, Shen et al. (2012) explored the possibility of deriving wind speed for longer ranges from the FSP of phased-array beam-forming HF WEllen Radar (WERA) systems. This work promoted the application of HF radar FSP in remote sensing of wave and wind parameters. In 2015, the FSP of a crossed-loop HF radar (OSMAR-S) was used to derive wave significant height by Zhou and Wen (2015). Later, Zhou et al. (2017) used the maximum FSP to estimate wind speed in a specific range, and obtained a result as good as that from the second-order inversion method. At the same time, Kirincich (2016) showed the wind field results of a broad-beam HF radar (SeaSonde) by the direct calibrated first-order scattering with the aid of an autonomous surface vehicle.

In this paper, we further analyze the effect of wind on the FSP. The analysis is based on HF radar scattering with a frequency of 13 MHz which means that the Bragg wave has a wavelength of 11 m. As the corresponding Bragg wave has a short wavelength, the Bragg wave will respond quickly to the local wind. Based on two semi-empirical models and a theoretical propagation attenuation model, we propose a space recursion method to extract the real-time wind field over the coastal ocean. The method can be applied to both broad-beam and phased-array systems and no modifications are required to adapt to either or both types of HF radar. This paper is structured as follows. Section 2 describes the theory and method. In Section 3 the field experiment and the analysis of the HF radar data are described. The theoretical and semi-empirical models are also formulated in this section. Model tests and wind vector mapping results are displayed in Section 4. Section 5 is the conclusion.

2. Theory and methods

The FSP of HF radar is influenced by many factors. Among these, the first-order cross section and propagation attenuation are the important parts related to the wind. We can write as

where $ {P}^{\left(1\right)} $ is the received FSP, $ u $ is the wind speed, $ \varphi $ is the wind direction, $ {\sigma }^{\left(1\right)}\left(\varphi,u\right) $ is the first-order cross section, $ A\left(\varphi,u\right) $ represents the absolute value of propagation attenuation and $ E $ is a constant related to other factors such as antenna and radar system. $ m $ and $ n $ are the range and angle cell number in a range-azimuth cell $ C(n,m) $ as shown in Fig. 1.

Figure  1.  Radar division of distance and angle.

DownLoad: Full-Size Img PowerPoint

The first-order cross section in cell $ C(n,m) $ is defined as

where $ {\sigma }^{\left(1\right)} $ is the first-order cross section, $ \omega $ is the Doppler frequency, $ {\omega }_{{\rm{B}}} $ is the first-order Bragg frequency, $ i=\pm 1 $ denotes the sign of the Bragg waves toward and away from the radar, $ {{\tilde k}}_{0} $ is the radar wave vector, $ {k}_{0}\!=\!\left|{{\tilde k}}_{0}\right| $ and $ \text{δ} (\cdot) $ is the Dirac delta function. $ S(\cdot) $ is the directional ocean wave spectrum and has the form:

where $ 2{k}_{0} $ is the scalar Bragg wave number, $ F\left(\cdot \right) $ represents the non-directional wave spectrum and $ G(\cdot) $ represents the wind spreading function (WSF). The WSF describes the directional energy propagation pattern of the ocean wave. When $ G\left(\varphi,u,2{{\tilde k}}_{0}\right) $ and $ A\left(\varphi,u\right) $ are obtained, the wind speed can be estimated by the $ F\left(\cdot \right) $. The FSP is assumed non-zero only when the sea is fully developed at a wind speed above a low wind cutoff: $ {U}^{2}\!\!>\!\!g/{k}_{0} $, where $ g $ is the gravitational acceleration and $ {k}_{0} $ is the radar wave number. Additionally, when the wind speed increases to a certain value, the Bragg wave associated with a given wave number tends to be saturated, which means the Bragg wave will no longer develop. The saturation property limits the range of wind speed estimation. For the 13 MHz HF radar, the cutoff wind speed is about $ 4 $ m/s (lower limit) and the saturated wind speed is about 13 m/s (upper limit).

2.1 Wind spreading function (WSF)

The WSF determines the accuracy of the wind direction inversion and is considered to be related to wind speed. So far, many WSF models have been proposed (Apel, 1994), and there is no broad consensus as to which WSF is optimal for a given situation. Among these WSFs, the cosine model and the sech model are applied widely by HF radar. For the cosine model, the focus is on the selection of the $ s $ parameter and the relationship between the $ s $ parameter and wind speed (Wyatt et al., 2006). Here, we used the WSF based on the cosine distribution as

where $ \theta $ is the relative angle between wind direction and the beam direction of arrival (DOA), and $ s $ is a parameter which is considered to be related to wind speed. Figure 2 shows the distribution of the WSF with different $ s $. The detailed relationship between parameter $ s $ and wind speed for the 13 MHz Bragg wave will be obtained from the radar data.

Figure  2.  Attenuation of first-order spectrum power (FSP) result from the wind spreading function (WSF) with different $s$.

DownLoad: Full-Size Img PowerPoint

According to the WSF, wind direction can be determined by the ratios of first-order Bragg peak intensities as

where $ {\omega }_{{\rm{B}}} $ represents the positive Bragg wave, $ -{\omega }_{{\rm{B}}} $ represents the negative Bragg wave and $ {{{R}}_{{{{\rm{Bragg}}}}}}$ represents the Bragg ratio. The wind direction estimation requires not only the Bragg ratio, but also the value of the parameter $ s $. Hence, the relationship between $ s $ and wind speed (the WSF model) needs to be established before the wind direction can be determined. The wind speed in cell $ C(n,m) $ cannot be obtained while estimating the wind direction and so the wind speed in the previous range cell is used as an approximation. Based on the wind direction and the WSF model, the FSP attenuation by the WSF can be calculated. In general, the FSP attenuation caused by WSF is large, especially when the DOA of the positive Bragg peak is from the down-wind direction as Fig. 2 shows.

2.2 Propagation attenuation model

HF ground-wave propagation losses across the ocean have been estimated for several decades. Barrick (1971) derived an expression for the effective impedance of a slightly rough, finitely conducting planar surface at grazing incidence for vertical polarization. He then gave an estimate of the propagation loss across the ocean using an empirical ocean height spectrum model. Forget et al. (1982) contrasted theoretical attenuation with actual HF radar attenuation by a monostatic experiment to evaluate the relative ground wave attenuation over the sea surface. His work shows that the radar relative attenuation has a good agreement with the attenuation of Barrick’s theoretical calculation. The propagation attenuation is closely related to the effective impedance of the ocean surface. The effective impedance at grazing incidence for vertical polarization can be written as

where $ {k}_{0} $ is the radar wave number, $ \theta $ is the angle between the DOA and wind direction, $ p $ and $ q $ are the radian wave number along the x and y directions. $ \varDelta$ is the normalized impedance defined as $\varDelta \!\!=\!\! {Z_s}/{Z_0}$, $ {Z}_{s} $ is the impedance of the surface sea in ohms and ${Z_0} \!\!\approx \!\!120\pi \varOmega$ is the impedance of free-space. $ W\left(p,q\right) $ represents a Neumann-Pierson spectrum, where $ C $ is a constant empirically estimated to be 3.05 m²/s⁵. Figure 3 shows the normalized impedance of the 13 MHz radar wave at different wind speeds. The normalized effective impedance increases with the growth of wind speed, it is larger in radial wind than cross wind. The higher impedance means greater attenuation for the same distance and radar frequency. Based on the normalized effective impedance, the attenuation can be calculated as Hill and Wait (1980). Figure 4 shows the attenuation under different sea states. We can see that the added attenuation caused by the rough sea surface is reduced almost linearly with distance within 50 km for both radial and cross winds. For the high SNR required for the analysis, we are only concerned with attenuation within 50 km. By fitting the attenuation curve by simple first-order linear fitting ($y\!=\!p_1x+p_{2}$), we obtain the requisite parameters in different sea states. The fitting parameters under different winds are plotted in Fig. 5. The parameters show the second order line distribution with the wind speed. By fitting as $p\!=\!a{u}^{2}\!+\!bu\!+\!c$, we obtain the parameters as $p_{1}\!=\!-0.000\;83{u}^{2}\!+\!0.001\;1u\!+\!0.103$ and $p_{2}\!\!= $$ \!\!-0.008\;2{u}^{2}\!+ \!0.010\;13u\!-\!0.064$, where $ u $ represents the wind speed in m/s. The effects of different beam angles as a function of wind direction are fitted by a surface which is presented as addloss=–0.57+0.02x+0.17u–0.004 6xu–0.010 3u², where $ x $ is the distance in km and $ u $ is wind speed. The relationship determines the added attenuation caused by the wind-direction related roughness. It should be noted that the propagation attenuation is calculated without consideration of the existence of dry or moist land between the radar site and the ocean cell. Land-sea interface can cause significant propagation loss and the loss has an angular dependence (Adams et al., 1984; Apaydin and Sevgi, 2010).

Figure  3.  Normalized surface impedance for 13 MHz radar wave under different wind speed, $ {0}^{\circ } $ represents radial wind and $ {90}^{\circ } $ represents cross wind. The solid and broken lines represent the real and imaginary parts of the impedance, respectively.

DownLoad: Full-Size Img PowerPoint

Figure  4.  Added loss for a rough sea under different winds. a. Under radial wind; b. under cross wind; c. under 20 kn wind speed; and d. under 30 kn wind speed.

DownLoad: Full-Size Img PowerPoint

Figure  5.  Fitting of the parameters under different radial winds (a), and discrepancy loss between radial wind and cross wind under different wind speeds (b).

DownLoad: Full-Size Img PowerPoint

For a monostatic HF radar, the received power can be written as

where $ \alpha $ is a constant term that includes the transmitted power, antenna gains, and system losses; $ \lambda $ is the electromagnetic wavelength; $ R $ is the distance of observation of the sea cell whose area is $ {\varDelta _s}$; $ \sigma $ is the radar cross section (RCS); and $ {F}^{2} $ is the attenuation factor for the energy relative to a one-way path between the radar and the sea cell. For HF radar, the $ {\varDelta _s}$ is proportional to the distance $ R $ since $ {\varDelta _s}$ is determined by the aperture angle of the radar beam and the distance $ R $, and hence the backscattered power is proportional to $ (1\!/\!{R)}^{3} $ (Forget et al., 1982). In addition, the amplitude of the receiving signal is proportional to the target distance due to time gate modulation in the radar transmitting frequency modulated interrupted continuous wave (FMICW) of the crossed-loop HF radar (Khan and Mitchell, 1991) as Fig. 6a shows. Figure 6 shows that the transmitted and received waveform are gated by a gating sequence with a $ 50 $% duty cycle. The modulation causes the magnitude of received signal to increase linearly with distance $ R $. Considering the scattering area and the influence of radar time gate control, the final power of the receiving signal is proportional to $ 1\!/\!R $ and the relationship between the basic free-space attenuation and the distance is shown in Fig. 6c. The final propagation losses are the basic free-space loss added to the sea surface loss as shown in Fig. 7.

Figure  6.  Modulation process of a gate signal to a receiving signal (a), effect of the gate signal on the radar echo power (b), and basic transmission loss for a monostatic 13 MHz OSMAR radar (c). Each range cell is $ 2.5 $ km. In a, A_T is the transmit pulse, G_R is the receive pulse, t₀ is the delay corresponding to distance R, A_R is the effective receive pulse for distance R, and T_p is the sweep period.

DownLoad: Full-Size Img PowerPoint

Figure  7.  Total transmission loss for a rough sea for the monostatic 13 MHz OSMSAR radar. a. Under radial wind; and b. under cross wind.

DownLoad: Full-Size Img PowerPoint

2.3 Block diagram of the wind field estimation

From the above equations we can see that $ {\sigma }^{\left(1\right)} $ and $ A $ are both related to the wind speed and wind direction. Hence, we propose a space recursion method to get the wind vector by the FSP. The overall procedure of the space-recursion method is illustrated in Fig. 8. The left side of the flow chart introduces the establishment of the FSP and WSF models, and the right side shows the extraction of wind direction and wind speed in each patch. In the establishment of the FSP model, the buoy wind speed can be replaced by the wind speed extracted from the second-order method. The wind-speed-related WSF model is a parameter $ s $ model for the cosine distribution established to estimate the wind more accurately than the $ s $ fixed cosine WSF. These models are important and their fitting accuracies determine the final wind estimation errors. In Section 3, we will describe the establishment of these models in detail. The wind speed field estimation by the FSP is the key to this method. The details of this method have these steps.

Figure  8.  Block diagram of the wind field estimation by HF radar.

DownLoad: Full-Size Img PowerPoint

(1) Build the second-order wind speed model, the WSF model and the FSP model as shown on the left chart of Fig. 8. Before applying the space recursion method, we need to establish the relationship between the maximum FSP and wind speed, namely, the FSP model. The FSP model is obtained by measuring the maximum FSP at the initial range cell which is about 5 km from the radar site. The short distance makes the effect of wind on propagation attenuation of the FSP negligible. The subsequent propagation attenuation is calculated from the initial range cell. In order to obtain more accurate wind direction attenuation and propagation attenuation, we also build a WSF model based on the cosine distribution and a theoretical propagation attenuation model. Both of these models are associated with wind speed. A month’s radar and buoy data are used to establish the second-order wind speed model, the FSP model and the WSF model between radar signal and wind speed.

(2) Obtain the approximate reference wind speed and the initial wind speed. The reference wind speed is used to judge the applicability of the method. When the reference wind speed is in the range of 4–13 m/s, the method can be applied in each patch or range-azimuth cell. The initial wind speed is used to estimate wind direction at initial range cell. The reference wind speed also can be used as the initial wind speed when it is obtained at the initial range cell. If there is no buoy wind speed at the initial range, the approximate reference wind speed and the initial wind speed $ {u}_{0} $ is determined by the second-order wind speed model as ${u}_{0}={f}_{2}\left({\int }_{\!\!-\infty }^{\infty }{\sigma }^{\left(2\right)}\left(w\right){\rm{d}}w\right)$, where $ {f}_{2} $ represents the second-order wind speed model. If there is a buoy wind speed measurement at the initial range cell, the buoy wind speed can be applied as the reference and initial wind speed $ {u}_{0} $ as shown in the gray square on the upper left of Fig. 8.

(3) Estimate the wind direction and wind speed at the initial range cell. At the initial range cell (in here we set the second range cell as the initial range cell), the propagation attenuation can be ignored. The wind direction is determined first in each patch. The wind direction and speed are obtained as

where $ {\varphi }_{0\left(m\right)} $ is the direction of radar echoes with angle number m, $ {\varphi }_{\left(2,m\right)} $ is the wind direction in the initial range cell $ C(2,m) $ with range number $ 2 $ and angle number $ m $ as shown in Fig. 1; $ {u}_{\left(2,m\right)} $ is the wind speed in cell $ C(2,m) $; $ {R}_{(2,m)} $ is the Bragg ratio in cell $ C(2,m) $; $ G $ represents the relative attenuation in decibels caused by the WSF, $ {p}_{\left(2,m\right)}^{+} $ is the positive radar FSP in cell $ C(2,m) $; $ {f}_{1} $ represents the maximum FSP model and $ s $ indicates the parameter $ s $ model for the cosine WSF model. The process is similar to the right part of Fig. 8, but the wind speed in the previous range cell is replaced by the initial wind speed, and the compensation of propagation attenuation is ignored.

(4) Estimate the wind vector in a non-initial range patch as shown in the right part of Fig. 8. Assuming that the wind speed for the WSF and propagation attenuation calculation is consistent with that of the previous range cell, we use it to compensate for the energy attenuation of the FSP due to the WSF and propagation. Then, the compensated FSP is applied to extract the final wind speed in this patch by the FSP model. This process is repeated across the radar scan area to get the wind vector in each cell by

where $ {\varphi }_{\left(n,m\right)} $ is the wind direction in cell $ C(n,m) $; $ {u}_{\left(n,m\right)} $ is the wind speed in cell $ C(n,m) $; $ {R}_{(n,m)} $ is the Bragg ratio in cell $ C(n,m) $; $ {p}_{\left(n,m\right)}^{+} $ is the positive radar FSP in cell $ C(n,m) $; and $ {A}_{\left(n,m\right)} $ is the propagation attenuation from the cell $ C(2,m) $ to the cell $ C(n,m) $.

3. Experiment and model parameters

3.1 Experiment and radar data

During February and March $ 2013 $, an HF radar named Ocean State Measuring and Analyzing Radar, type S (OSMAR-S) was installed on the west coast of the Taiwan Strait to monitor the sea state. The geographic map is illustrated in Fig. 9. The location of the OSMAR-S was $ {23.749}^{\circ } $N, $ {117.597}^{\circ } $E. The OSMAR-S is a compact, broad-beam all-digital HF radar system. It adopts a linear FMICW waveform with a center frequency of 13 MHz and a bandwidth of 60 kHz. The range resolution is $ 2.5 $ km correspondingly. The specific radar parameters used for this test are listed in Table 1. The reliability and accuracy in current and wave measurement of the OSMAR-S have been validated in Lai et al., (2017a, b). Meanwhile, three in-situ wave-rider buoys, namely, A ($ {23.783}^{\circ } $N, $ {118.033}^{\circ } $E), C ($ {23.67}^{\circ } $N, $ {117.67}^{\circ } $E) and E ($ {23.417}^{\circ } $N, $ {117.917}^{\circ } $E) acquired wind data every half-hour. The distances between the radar station and the three buoys were 44 km, $ 10.5 $ km, and $ 48 $ km, respectively. During the experiment, the wind predominantly blew from the northeast with speeds that varied between 0 m/s and 16 m/s as shown in Fig. 10. Buoy E had a relatively high wind speed and a relatively stable wind direction, while Buoy C had a relatively low wind speed. In general, the experimental conditions are ideal because the winds are stable and the effects of wind fetch (infinite wind fetch for northeasterly wind) and water depth (above half of the Bragg wave length as shown in Fig. 9) are negligible.

Figure  9.  Map of the HF radar deployed on the Fujian coast of the Taiwan Strait. Radar station is marked by the star named XIAN, and the three buoys are marked by dots.

DownLoad: Full-Size Img PowerPoint

Table  1.  Radar parameters of the OSMAR-S

Parameters	Value and setting
Center frequency/MHz	13
Bandwidth/kHz	60
Transmit antenna	monopole
Receive antenna	cross-loop/monopole
Sweep period/s	0.38
Average power/W	100
Range resolution/km	2.5
Coherent integration time/min	6.5
Normal direction/(°)	100
Transmitted waveform	FMICW pulses
Technique of azimuthal resolution	direction finding

 | Show Table

DownLoad: CSV

Figure  10.  Wind rose figure of Buoys A (a), C (b) and E (c).

DownLoad: Full-Size Img PowerPoint

We built the models based on the radar data collected from February 1 to 28, 2013. To ensure good quality, data with SNR$ < $5 dB are removed as are any data with radial current velocity above the expected maximum (1 m/s) for this region (Lai et al., 2017). The spectral region on the side with a stronger Bragg peak is chosen to estimate the DOA by Multiple Signal Classification (MUSIC) direction-finding algorithm (Schmidt, 1986) using the ideal antenna pattern. Figure 11 shows the times of positive and negative first-order peaks in each cell for 6 202 frames of data. From Fig. 11, we can see that most DOAs are near the antenna normal direction (100° from the north). The number of viable positive first-order peaks is more than twice that of the negative ones because the values of the positive FSP are always larger than those of the negative ones. This means that the negative peaks from the southwest are lower in amplitude than the positive peaks from the northwest of each range annulus. To avoid the power difference caused by the distortion of antenna pattern in different directions, the FSP model is only based on the positive FSPs of the monopole-antenna.

Figure  11.  Number of DOA solutions in each measurement bin, where the total number of measurements is 6 202. a. For positive FSPs; and b. for negative FSPs.

DownLoad: Full-Size Img PowerPoint

3.2 Propagation attenuation analysis of radar data

In this section, the actual HF radar propagation attenuation is calculated from one-month radar data. Figure 12a shows the maximum FSP in different wind speeds and ranges. The maximum FSP is assumed to come from the upwind direction and the effect of WSF can be ignored. From Fig. 12a, we can see that the radar echoes seem to reach the maximum when the wind speed is about 10 m/s. Subsequently, the maximum FSP decreases with the increased wind speed, especially in the far range cell. It results from the increased sea surface roughness. Figure 12b shows the maximum FSP loss (relative to the maximum FSP in range cell 2) under different wind speeds. The maximum FSP loss is the total propagation attenuation and seems to be attenuated linearly with the distance when the wind speed is between 3 m/s and 12 m/s. With the increase of the range cell, the differences between the attenuation under different wind speeds increase gradually. In range cell 23, the attenuation difference between the maximum and the minimum wind speed is up to about 13 dB, which is consistent with the theoretical calculation. Figure 13 shows both the theoretical and the radar extracted propagation attenuation. Basically, the difference is less than 5 dB. We apply the theoretical propagation attenuation model to compensate for the propagation attenuation.

Figure  12.  The maximum FSPs versus wind speeds for different range cells (a); and the maximum FSP loss (relative to the maximum FSP in range cell $ 2 $) versus range cells under different wind speeds (b). Each range cell is $ 2.5 $ km.

DownLoad: Full-Size Img PowerPoint

Figure  13.  Comparison of the theoretical and the radar extracted propagation attenuation relative to the range cell 2 with different radial wind.

DownLoad: Full-Size Img PowerPoint

3.3 Wind spreading function model

The principle of HF radar wind direction inversion is based on the first-order Bragg ratios as Eq. (5) shows. The WSF is dependent on wind speed through the parameter $ s $. Its compensation is an important part of the FSP attenuation compensation. Hence, we extracted a WSF model by establishing a parameter $ s $ model based on one month’s radar data with the aid of buoy wind data by

where $ s $ is the parameter of cosine model as shown in Eq. (4), $ {R}_{{\rm{B}}} $ is the Bragg ratio, and $ {\theta }_{0} $, $ {\theta }_{\omega } $ are the DOA and the buoy wind direction, respectively. In Fig. 14, the radar parameter $ s $ is calculated under different buoy wind speeds. When the wind speed is less than $ 13 $ m/s, the parameter $ s $ increases with the growth of wind speed and it decreases when wind speed is greater than $ 13 $ m/s. Compared with $ s $ from the literature (Hasselmann et al., 1980), the variation of radar extracted $ s $ with wind speed is flatter as a function of wind speed. According to this distribution, a second-order linear fitting is used to fit this distribution, and we get the relationship that $s\!=\!-0.010\;6{u}^{2}\!+\!0.256\;4 $$ u\!+\!1.884\;5$, where u is the wind speed. The correlation coefficient (R) and root mean square error (RMSE) are 0.76 and 0.21, respectively. As mentioned above, the wind speed in the WSF model is approximated by the wind speed in the previous range cell.

Figure  14.  Parameter s versus wind speed for the maximum first-order peaks.

DownLoad: Full-Size Img PowerPoint

3.4 Second-order wind speed model

We apply the second-order scattering in the initial range cell to obtain the preliminary reference wind speed. The preliminary wind speed can be used to determine whether the first-order peak is saturated and to initialize the wind vector mapping process when the buoy data are missing. Figure 15a shows the time series of C buoy wind speed and the second-order spectrum integration at range cell 2. These two variables show a strong correlation with a R of $ 0.85 $. We plot it in a scatter plot and fit it by a first-order linear fitting which has an RMSE of 2.12 m/s and an R of 0.88. The expression of the linear fitting is $ u\!=\!0.003\;7\!\times\!P_s- $$ 27.57 $ where $ P_{\rm{s}} $ is the second-order spectrum integration and u is the wind speed. According to the Pierson-Moskowitz (PM) non-directional wave spectrum (Pierson and Moskowitz, 1964) and the radar echo spectrum (Barrick, 1977), a first-order linear relationship between wind speed and second-order spectrum integration can also be developed:

Figure  15.  Wind speed and the second-order integration versus time (a); scatter plot and linear fitting result (b).

DownLoad: Full-Size Img PowerPoint

where $ H_{{\rm{s}}} $ is the significant wave height and $ W\left(\omega \right) $ is a weighting function. Although this relationship is an approximate and depends on a number of assumptions, it demonstrates the reason for the linear relationship between the second-order integration and the wind speed. It can be seen that the correlation between the two is better at high sea state and worse at low sea state for the absence of second-order spectrum.

3.5 First-order spectrum power (FSP) model

The FSP model is the key part for wind speed estimation. After the energy compensation for the FSP, the final wind speed can be extracted directly by the FSP model. To avoid the uncertainty caused by WSF, we track the maximum FSP of the monopole at the initial range cell. Because the model is based on the radar echoes at the initial range cell (about $ 5 $ km from the radar station), the added attenuation of FSP from the rough sea under different wind speeds can be ignored as Fig. 4a shows. According to Eq. (2), the Bragg spectral value is a sample of the directional wave height spectrum at the Bragg frequency. The directional wave height spectrum can be modeled as the product of the non-directional wave height spectrum and the WSF as Eq. (3) illustrates, which means that the wind speed is closely related to the maximum FSP when the wind speed is below 13 m/s. Figure 16a shows the maximum FSP and the buoy wind speed in time series. To reduce the effect of noise, a five-point smoothing filter is applied to the time series of the maximum FSP. These two variables have a strong correlation with an R of 0.88. Their scatter plot shows a power distribution, and it fits the model as

Figure  16.  Wind speed and the maximum first-order peak versus time (a); and scatter plot and fitting result (b).

DownLoad: Full-Size Img PowerPoint

where $P_{\max} $ is the maximum FSP in dB and $ a $, $ b $ and $ c $ are constant parameters that are determined by the training data. The a, b and c are $ -1.096\!\times\!{10}^{7} $, 29 and 119, respectively. The maximum FSP model can provide a good fit to the data with an R of 0.9 and an RMSE of 1.54 dB as Fig. 16b shows. The fitting results of the model are better when the wind speed is in the range of 4 m/s to 13 m/s. The FSP and wind speed are more relevant at low sea state and irrelevant at high sea state for the saturation of the FSP.

4. Model test and mapping

To verify the reliability of the model, we test the three models with the radar and buoy data in March. The propagation attenuation model is consistent with the theoretical model and can be considered to have high reliability, hence we no longer test it.

4.1 Second-order wind speed model test

When there is no buoy data, the second-order model based wind speed can be used as an approximate reference and an initial wind speed. Figure 17 shows the wind speed comparison between the buoy and that obtained by the inversion of the second-order radar backscatter. The R between the two is 0.68. Around March 17, the radar radio interference was high, which we believe resulted in the overly-large wind-speed estimation. Except for this period of interference, the second-order integral model appears to perform well throughout March.

Figure  17.  Buoy wind speed versus the second-order integration (a); and buoy wind speed versus the radar wind speed from the second-order model (b).

DownLoad: Full-Size Img PowerPoint

4.2 First-order FSP model test

Accuracy of the FSP model is very important for the reliability of the method. Therefore, it is necessary to test the FSP model. Figure 18 shows the comparison between the buoy wind speed and the radar wind speed inverted from the FSP model. To reduce the effect of noise, a 17-point smoothing filter is applied. The R between these two is 0.80 with an RMSE of 2.84 m/s. On March 7 and 8, ships and interference are in the radar first-order spectrum area, which leads to the higher errors in wind speed estimation. The results demonstrate the applicability of the FSP model.

Figure  18.  Comparison of the wind speed and the radar wind speed from the first-order model.

DownLoad: Full-Size Img PowerPoint

4.3 Wind direction test

Figure 19 shows the radar inverted and buoy wind direction at A for wind speed above 6.5 m/s. When wind speed is above 6.5 m/s, the wind directions are more stable and reliable. Figure 14 shows the WSF model used in the wind direction inversion. To test the reliability of the model, we use the buoy wind speed in the formula and a prior reference wind direction is used to remove the directional ambiguity. In this case, the radar tends to underestimate the wind direction and the bias of the wind direction estimation is $ {22.83}^{\circ } $. It should be noted that radar data in March are more frequently disturbed by ship and radio interference than the data in February. The interference makes the error rate of wind direction high.

Figure  19.  Sample comparisons of the in-situ measured and estimated wind direction at Buoy A in March. Note the sample numbers do not represent a time series but are individual, independent estimates of half an hour time period with wind speed above $ 6.5 $ m/s.

DownLoad: Full-Size Img PowerPoint

4.4 Mapping result

Based on the three models we can extract the wind field (Fig. 20) from the first-order spectrum by the space recursion method as shown in the flow chart above. The wind speed extracted from the second-order spectrum is used to initialize the wind vector mapping process. To obtain a stable wind map, a Gaussian smoothing filter is applied, which is reasonable and widely used due to the typical continuity of the wind field. The estimated wind is predominately blowing from northeast to southwest along the strait, which is in line with the in-situ measurements. The differences between the estimated wind speeds are small in most areas. The differences between the estimated wind speeds are small in most areas. Since the DOAs of the positive first-order peaks are more from the normal direction as shown in Fig. 11, the wind estimates are more satisfactory in this region. Whereas, the wind estimations in the north of Buoy A are different from those in other parts of the wind field. The estimated wind directions in the north area seem perpendicular to the radar beam and they show a circumferential pattern. In this area, the Bragg ratios are small and tend to estimate wind direction to be perpendicular to the radar beam. The small Bragg ratios may be caused by wave refraction in the near-shore area, or the uncertainties in the actual distribution of the wave power along the direction, or the error of estimated DOAs. Because there is no additional wave measurement, it is difficult to separate small Bragg ratios from errors. For the wind speed, most estimations are consistent with the buoy’s wind speed, especially near the antenna normal direction area. Due to the island-caused additional propagation loss, the estimated wind speed at A and the north area of A is small. The radar-derived wind speed at long distance (above 18 range cell) is also relatively small because the compensation of the propagation attenuation is not enough. The error in propagation attenuation model has little effect on the near distance element (below range cell 10) as shown in Fig. 13. From Fig. 13, we can see that for the propagation attenuation within range cell 10, the difference between the theoretical and the actual attenuation is less than $ 5 $ dB for different sea states. With the increase of the distance, the wind speed error accumulates gradually and a more accurate propagation attenuation model is needed. The radar-derived wind speeds at the west side of Buoy E are relatively large because of the fewer DOAs and the overcompensation of the WSF.

Figure  20.  Final wind field inversion results at 7:30, March 4, 2013.

DownLoad: Full-Size Img PowerPoint

4.5 Discussion

By the FSP model, the radar estimated wind speed has an R of 0.8 and an RMSE of 2.84 m/s with the buoy. The wind direction estimated from the WSF model has an RMSE of $ {22.83}^{\circ } $ at wind speeds greater than 6.5 m/s. The error rates of the wind inversion are generally comparable to the results of recent wind extraction studies. Similar satellite-based wind extractions to buoy wind observation achieved the RMSE $ 1.7\!-\!2 $ m/s and $ {40}^{\circ }\!-\!{50}^{\circ } $ for speed and direction, respectively (Carvalho et al., 2014; Zhang et al., 2016). For HF radar, Wyatt et al. (2006) presented 6 month HF radar wind speeds and wind directions with RMSE of 5.188 m/s and $ {23.13}^{\circ } $ compared with the model products between a big range of wind directions. Green et al. (2009) used the position shift of second-order peaks to estimate the wind speed with error rates about $ 20 $%. Kirincich (2016) gave wind direction estimate with RMSE between $ {35}^{\circ } $ and $ {40}^{\circ } $ and wind speed estimate with RMSE ranging from 1.4 m/s to 2.7 m/s. It is important to note that only Kirincich’s method can provide wind field mapping for a broad-beam HF radar.

For the space recursion method, the propagation attenuation calculation, the WSF and the FSP models are the three main parts. The second-order integral model is only applied at the initial range cell to get a rough sea-state estimation, which is used to judge whether this method can be applied or to initialize the wind speed inversion. For a 13 MHz HF radar, the method can be applied when the wind speed is in the range of 4–13 m/s. When the wind speed is greater than 13 m/s, the change of the first-order peak with the wind speed is not obvious, and the accuracy of estimated wind speed will decrease.

The compensation of propagation attenuation is the foundation of using the FSP. Here, we have considered the influence of wind-driven sea surface roughness on the echoes’ propagation attenuation. Therefore, wind information is required to calculate the total propagation attenuation. In practical application, we use the wind of the previous range cell as the approximate value to estimate the influence of wind on echo attenuation. Based on the theoretical calculation of the $ 13 $ MHz transmission loss for a rough sea, we built a simple wind-related propagation attenuation model for the compensation. The model’s reliability is verified within a distance of 50 km. At longer distances, the model needs additional adjustments. The effect of dry or moist land on the model is neglected. The propagation loss over a land is about twice that over a sea (Adams et al., 1984) and the land-sea interface will also result in dramatic propagation loss (Apaydin and Sevgi, 2010). Here, there are some islands on the northeast of the radar site. These islands will decrease the FSP in the northeasterly direction. The decreasing of the FSP leads to a small wind speed estimation as shown on the northeast of Fig. 20.

The WSF model not only has great influence on the accuracy of wind direction estimation but also is the key part for wind speed compensation. As Fig. 16b shows, in the range of wind speed from $ 1 $ m/s to $ 15 $ m/s, the average first-order power only varies by about $ 15 $ dB. However, for the WSF, the angle changes of $ {180}^{\circ } $ can produce nearly $ 35 $ dB of energy attenuation as shown in Fig. 2. Therefore, the WSF and Bragg ratio are very important for the wind estimation. When the radar echoes are from the upwind direction (i.e., $ \theta \!=\!0 $° according to Fig. 2), the wind direction estimation error will not generate large attenuation of the FSP, which means a small wind speed estimation error from the WSF. At present, the wind direction estimation error of HF radar is roughly $ {20}^{\circ }\!-\!{40}^{\circ } $, which will cause less than 5 dB attenuation around the upwind direction according to the WSF as Fig. 2 shows with $ s\!=\!3 $.

The WSF model can be replaced by the existing oceanographic model, but the FSP model must be trained because there is no theoretical model. Moreover, the FSP model is closely related to the radar system and to the impedance of the sea surface (which can be affected by the land-sea interface or the seawater salinity), which means that for a different radar system or in a different sea area that the model needs to be trained again. Assuming deep water of the radar location, the training buoy is as close as possible to the radar. As shown in Fig. 7, the closer the buoy is to the radar, the less the difference between the propagation attenuation under different sea states and the more accurate the FSP model. Generally, the buoy location is preferably within 10 km of the radar.

In the process of using the method, we must consider the saturated wind speed range of the corresponding Bragg wave, especially in extreme cases like coastal storms. In these cases, the wind speed easily exceeds the range of saturated wind speed and a radar with lower frequency is needed. It is also important to note that the FSP is susceptible to radio frequency interference, background noise, etc. Hence we should judge these factors and eliminate their influence before using the method. When the method is applied to a broad-beam radar, we should make sure that the DOAs are estimated correctly. Multiple DOA condition may occur when currents are parallel to the radar beam. When some wind conditions are not seen during the training period, the method is still usable as long as the fitting can convergence. However, loss of large amounts of training data will affect the accuracy of the model and hence the final result.

5. Conclusions

HF radar has been an effective tool in sea state remote sensing. However, wind mapping with HF radar is still a problem, especially for a broad-beam HF radar. In this paper, we give a useful space recursion method to map the wind vector by first-order Bragg scattering. With the help of one-month of buoy data, the relationship between the FSP and the wind has been fully analyzed. The FSP and WSF models are established in the near shore for wind mapping. The models are tested using one month radar data and their performances are supported by the test results. The wind speed estimated from the FSP model and that from buoy have an R of $ 0.8 $ and an RMSE of $ 2.84 $ m/s. The wind direction estimated from the WSF model has an RMSE of $ {22.83}^{\circ } $ when wind speed is more than $ 6.5 $ m/s. The FSP and WSF model can be applied to other HF radar systems and in other sea areas after a calibration. The final wind vector results validate the effectiveness and practicability of the space recursion method for the broad-beam HF radar. Also, it can be used for narrow-beam, phased-array HF radar for wind vector mapping. The FSP space recursion method can cover a wider area than the second-order method when the wind speed is in the unsaturated range of the first-order Bragg spectrum. In practical applications, these two methods can be combined in wind field extraction. When the SNR of the second-order spectra is below a set threshold, the corresponding FSP and WSF model can be established by the wind speed extracted from the second-order echoes, and the FSP space recursion method can be applied in the outside areas. The amplitude of the first-order peak is the key parameter of this method and a smoothing of the FSP is necessary to ensure the performance due to the effect of noise. At the same time, the accurate direction finding or DOA estimation also can improve the result of the FSP based space recursion method. In the future, we will use more data and experiments to verify the reliability of this method.