Abstract

New signal processing algorithms for the Airport Surveillance Radar-9 (ASR-9) Weather Systems Processor (WSP) are introduced. The Moving Clutter Spectral Processing for Uneven-Sampled Data with Dealiasing (MCSPUDD) algorithm suite removes isolated moving clutter targets and corrects aliased velocity values on a per-range-gate basis. The spectral differencing technique is applied to the low- and high-beam data to produce a dual-beam velocity estimate that is more accurate than the current autocorrelation-lag-1-based approach. Comparisons with Terminal Doppler Weather Radar (TDWR) data show that estimate errors are reduced by 8%, 15%, and 15% for the low-, high-, and dual-beam velocities, respectively.

1. Introduction

Microburst-induced airline accidents in the 1970s and 1980s prompted the Federal Aviation Administration (FAA) to commission the development of wind shear detection systems. Resulting from this effort were the Low Level Wind Shear Alert System (LLWAS), a network of anemometers placed around airport runways and corresponding approach and departure paths (Wilson and Gramzow 1991), and the Terminal Doppler Weather Radar (TDWR), a stand-alone unit designed specifically for detecting hazardous weather phenomena in the terminal area (Michelson et al. 1990). The TDWR is a high-performance weather radar that is still considered the gold standard for terminal-area wind shear detection. There are 45 operational TDWRs in the United States. However, it was deemed too expensive to be deployed at airports that did not have a large rate of microburst occurrence or high air traffic density. A cost-effective alternative solution was proposed, which was to piggyback a signal processing system on the Airport Surveillance Radar-9 (ASR-9) (Taylor and Brunins 1985) that would provide wind shear detection capability at a subset of airports not served by TDWR or LLWAS. This system, the Weather Systems Processor (WSP) (Weber and Stone 1995), was ultimately deployed at 34 airports. Field evaluations have shown the WSP to adequately detect wind shear (Weber et al. 1996).

The difficulty of processing ASR-9 data for low-altitude wind shear observation is twofold. First, unlike a radar designed for weather surveillance like the TDWR or NEXRAD, the ASR-9 does not have a narrow antenna beam. Instead, it has a beam pattern that is specially shaped to sweep through the terminal airspace volume rapidly with the desired energy density on target. Thus, the pattern is relatively narrow in azimuth (1.4°) but broad in elevation (4.8°) with a cosecant-squared rolloff on the top side. Consequently, discrimination of the thin microburst outflow layer from ground clutter below and storm core aloft is harder. Second, because the antenna rotates much more quickly (75° s−1) than on a weather radar, the coherent processing interval (CPI) per azimuthal sector has to be short, which degrades the base data estimate quality and renders clutter filtering more problematic.

The current operational WSP signal processing algorithms mitigate these problems by making use of the low and high antenna beams available on the ASR-9, and by processing longer overlapped segments of data [extended coherent processing interval (ECPI)]. However, the processor and memory capacity of the original digital signal processor (DSP) limited the scope of techniques that could be implemented (Weber 2002).

With the first-generation WSP hardware nearing its end of designed life, we are currently developing a sustainable replacement for the future. As part of this technology refresh, the signal processor capability will be enhanced, which gives us an opportunity to likewise upgrade the processing algorithms. This paper reports on the new algorithms and presents data that quantitatively demonstrate improvement in velocity estimate quality over the legacy algorithms.

2. Brief review of legacy signal processing algorithms

The ASR-9 has a reflector antenna with two feed horns. The feeds have associated elevation patterns that overlap but are displaced in elevation angle. The low beam is always used for transmission, whereas both beams are used on reception. For aircraft target surveillance, the high beam receive channel is used at short range (typically to 28 km) in order to reduce ground clutter interference. The low-beam receive signal is used at farther ranges to maintain low-altitude coverage. The ASR-9 also has the capability to transmit either vertically polarized or right-hand circularly polarized signals. On receive, both copolarized and cross-polarized (horizontal in the linear case and left-hand sense in the circular case) signals are available. The linear mode is used during clear weather. When significant precipitation is detected in the coverage area, the polarization is switched to circular mode in order to reduce weather clutter in the target channel.

For optimal weather observation, the copolar received signal is needed during linearly polarized transmission and the cross-polar signal is needed during circularly polarized operation. Hydrometeors tend to be spherical and reflect most of the circularly polarized energy in the opposite sense. Therefore, the WSP taps into both the copolar and cross-polar radio frequency (RF) receive paths. However, since the WSP has only one receiver channel, its input is switched between the low and high beams after each “scan.” (A scan in this case was defined to be a 465.5° rotation in order to keep the switching transients from piling up at any particular azimuth angle.) This means that simultaneous coherent processing of low- and high-beam in-phase and quadrature (I&Q) time series data is not possible, ruling out techniques such as dual-beam interferometry. The I&Q data from the aircraft surveillance channel are available to the WSP, but they are generally not used because of the reduced data quality (e.g., lower dynamic range). The only exception is during linear polarization transmission when the far-range low-beam copolarized data are only available from the target channel. See Weber (2002) for further details on this and subsequent legacy signal processing topics.

The aircraft surveillance CPI consists of 10 pulses at a short pulse repetition interval (PRI) and then at least 8 pulses at a long PRI. If the end of the azimuth sector boundary (1.4° width) is not reached at that point, then fill pulses spaced at the long PRI are transmitted until the boundary is reached. The WSP ECPI adds to this basic CPI the long-PRI set of pulses from the previous azimuth sector. Thus, the ECPIs overlap each other across the azimuth sector boundaries by the length of the long-PRI pulse set. If there are fill pulses, then the ECPI is trimmed at the start so that it always consists of 27 pulses. The exact values of the PRIs can vary from site to site, but the ratio of long to short PRI is kept constant at 9:7.

Once the 27-pulse ECPI is collected at a given azimuth sector and range gate, and the amplitude is properly normalized for the sensitivity time control (STC) function, it is processed according to the block diagram in Fig. 1. First, adaptive clutter filtering is performed. The gist of this step is the comparison of data filtered at four clutter suppression levels (all pass, 20, 40, and 60 dB) with clutter residue maps (CRMs) at the same suppression levels collected during clear weather. The algorithm is designed to select the least attenuating filter that generates a weather-signal-to-ground clutter ratio that is big enough for quality velocity estimation. The clutter filters are of finite impulse response (FIR) type, yielding an excellent balance of magnitude response and phase linearity for block-staggered PRI pulse trains (Chornoboy 1993; Cho and Chornoboy 2005). If the filtered data signal power is not sufficiently higher than the corresponding clutter residue, then a “filter fail” flag is set to indicate that the weather-signal-to-clutter ratio is too small for wind shear detection processing. The adaptive clutter filter also contains an anomalous propagation (AP) detector. If the all-pass filter was chosen, but further clutter filtering attenuates the signal more than is expected for weather, then the range gate is flagged as contaminated by AP.

Fig. 1.

Legacy WSP signal processing block diagram. I and Q are in-phase and quadrature data, respectively; R0 and R1 are the autocorrelation lags 0 and 1, respectively; LZ is low-beam reflectivity; HZ is high-beam reflectivity; LV is low-beam velocity; HV is high-beam velocity; DZ is dual-beam reflectivity; DV is dual-beam velocity; LF is low-beam flag; and HF is high-beam flag.

Fig. 1.

Legacy WSP signal processing block diagram. I and Q are in-phase and quadrature data, respectively; R0 and R1 are the autocorrelation lags 0 and 1, respectively; LZ is low-beam reflectivity; HZ is high-beam reflectivity; LV is low-beam velocity; HV is high-beam velocity; DZ is dual-beam reflectivity; DV is dual-beam velocity; LF is low-beam flag; and HF is high-beam flag.

After clutter filtering, autocorrelation lags 0 and 1 (R0 and R1, respectively) are computed separately for the short- and long-PRI pulse sets. If the ratio of to deviates from unity by a certain threshold, then the larger R0 is deemed to be contaminated by out-of-trip signals and the appropriate flag is set (there is, however, a running average of this condition maintained across scans to avoid spurious flagging of out-of-trip contamination). If both PRI sets are free of range-aliased signals, then R0 and R1 are computed from the weighted average of and , respectively. Otherwise, R0 and R1 are assigned values from the clean PRI lags.

The lags are smoothed over time by a first-order recursive averaging filter across scans; that is, every time there are new lags data, the new value times a weight, A, is added to the current average value times (1 − A). The default value of A is 0.5. This lags buffer is maintained for each beam (low and high) at every azimuth and range position.

The reflectivities for the low and high beams, ZLo and ZHi, respectively, are then calculated from the R0s by subtracting the system noise, multiplying by the range squared and the ASR-9 radar constant, then converting to decibel units (Weber 2002). The low- and high-beam velocities are computed from the standard pulse-pair equation (Doviak and Zrnić 1993), except that the effective PRI (TEff) is calculated based on whether the short, long, or weighted average PRI sets were used to generate the R1s.

The dual-beam reflectivity is a weighted average of ZLo and ZHi, where the weights are range dependent. This dependency is a site-adjustable parameter, but typically it is set such that the high beam is wholly used at short range and the low beam is entirely used at long range.

The dual-beam velocity estimate was formulated in an attempt to leverage the overlap in the low- and high-beam patterns to sharpen the focus on the velocities at the bottom part of the low beam. It is the pulse-pair velocity estimator applied to a weighted linear difference of the low- and high-beam lag-1 autocorrelation estimates, R1Lo and R1Hi, respectively:

 
formula

where λ is the radar wavelength and W is a site-adaptable weighting parameter (Weber 2002). The symbol “∠” denotes taking the phase angle of the complex argument in the interval −π to π.

Finally, in the postprocessing step, censoring (based on the data quality flags) and 3 × 3 (range × azimuth) median filtering is performed on the dual-beam reflectivity and velocity fields.

3. Description of the new signal processing algorithms

We targeted two areas of data quality improvement with the new algorithm development. First, velocity dealiasing. The legacy processing algorithm makes no attempt to dealias velocity estimates, even though data from two sets of PRIs are available for this purpose. Consequently, velocity magnitudes greater than the Nyquist span of the long PRI (about 28 m s−1) are corrupted by aliasing. For reference, the requirement for the dealiased velocity interval is ±40 m s−1 for the TDWR and ±50 m s−1 for the NEXRAD. Furthermore, there is no provision for flagging and censoring aliased velocity data except for the relatively crude standard-deviation-based spatial filtering in the postprocessing phase. Since wind shear detection relies mainly on accurate and reliable radial velocity estimates, it is important to correct aliased velocities.

Second, the autocorrelation-based dual-beam velocity estimate relies on a somewhat arbitrary weighting parameter, W as seen in (1). In theory, there is an optimal W for every situation, but there is not enough information available to adaptively adjust W in real time (Weber and Noyes 1988). Also, there is no direct way to dealias the dual-beam velocity using this approach.

Our solution to these problems is to compute dealiased band-limited Doppler spectra for the low and high beams and then to employ a spectral differencing method to generate the dual-beam velocity. Figure 2 shows the outline of this new strategy. It follows the same processing procedure up to the clutter filter and autocorrelator and then deviates to perform velocity dealiasing and transformation to the Doppler spectral domain. More detail is given in the following subsections.

Fig. 2.

New WSP signal processing block diagram.

Fig. 2.

New WSP signal processing block diagram.

a. Velocity dealiasing

As described in the fifth paragraph of section 2, the short- and long-PRI sets in each ECPI are checked for out-of-trip contamination. If both sets are clean, then velocity dealiasing can be attempted. First, compute velocities corresponding to and :

 
formula

where TShort and TLong are the mean PRIs of the short- and long-PRI sets, respectively. The dealiasing can now be done via a clustering algorithm (Trunk and Brockett 1993) or a rule-based table lookup procedure (e.g., Torres et al. 2004). We opted for the latter due to its faster execution time. The steps are given below.

Generate a table that relates the theoretical velocity difference, VDiff = VShortVLong, to the aliasing (or folding) condition for VShort and VLong. For example, if VDiff = 0, then neither VShort nor VLong is folded, so and are both zero. If VDiff = −2, then = 0 and = −1, where = λ/(4TShort,Long) is the Nyquist velocity associated with each PRI set. The table is only generated to cover ±VMax, the desired dealiasing span (currently set to ±42 m s−1). It can be pregenerated offline, since it is only dependent on the mean PRI values and VMax.

Note that the Chinese remainder theorem (e.g., Ding et al. 1996) tells us that the ASR-9’s integral PRI ratio of 9:7 yields an unambiguously unfoldable span of ±9 = 7, which is about ±235 m s−1 (depending on each site’s radar wavelength and PRI values). By restricting the dealiasing operation to a more reasonable ±VMax span, we reduce the probability of incorrect dealiasing due to noisy velocity estimates.

Next, we find the theoretical value of VDiff in the lookup table that is closest to the measured VDiff. Use the corresponding folding values to generate the candidate dealiased velocity values for each PRI set:

 
formula
 
formula

However, if either of these candidate dealiased velocity magnitudes exceeds VMax, then find the next closest VDiff match in the lookup table. Repeat these steps if needed until both unfolded velocities are within ±VMax.

Now compute the dealiased velocity as a weighted average,

 
formula

where and are the number of pulses averaged to get VShort and VLong, respectively. Repeat the unfolding and dealiasing steps of (3)(5) for the next closest VDiff match in the lookup table. The corresponding dealiased velocity will be the “alternate” dealiased velocity, , which will be used in the false dealias correction step. If there were no more VDiff values in the table after the primary VDealias was computed, then mark this condition by placing a corresponding flag value in .

Finally, convert the dealiased velocity (and alternate dealiased velocity, if available) to the normalized mean frequencies,

 
formula
 
formula

where Va = λ/(4Tu) and Tu = TShort/7 = TLong/9 is a uniform subsampling interval. If one of the PRI sets was contaminated by out-of-trip signals and thus dealiasing was not possible, then assign the normalized mean frequency based on the “raw” velocity estimate from the clean PRI set. The normalized mean frequency will be used to phase shift the time series so that the signal is centered in the band-limited spectrum window within the ±Va span. But, first, we need to correct any falsely dealiased velocities.

b. False dealias correction

Because of the inherent variance in real weather radar data, false velocity dealiasing is unavoidable for some fraction of cases no matter what technique is used. The human eye can usually detect such errors, because of the available contextual information in space and/or time. Similarly, automated algorithms can also detect and correct such errors based on continuity. We adapted the following algorithm from a previously developed procedure for TDWR (Cho 2005).

First, make a copy of one radial’s worth of the dealiased normalized mean frequency data—call this sCp. At each range gate j, check for a false dealias condition if velocity dealising was performed and there was an alternate dealiased normalized frequency, sAlt(j), available. Otherwise, skip this gate.

Second, collect the neighboring sCp values from gates jNHalfWin to j − 1 and j + 1 to j + NHalfWin, but do not go below 1 or above the available number of range gates. Currently, NHalfWin is set to 2.

Third, eliminate any neighboring sCp values in this collection that were bad—that is, if the out-of-trip flag, clutter filter fail flag, or AP flag was set. If there are two or more neighboring sCp values left, then compute the median of those values, sMed(j). Otherwise, skip this gate.

Finally, if , then replace s(j) with sAlt(j).

c. Doppler spectral processing

To compute the Doppler spectrum of the WSP’s irregularly sampled data, we begin by following Chornoboy and Weber’s (1994) band-limited interpolation approach; that is, we represent the dual-PRI time series as an incomplete sampling of a hypothetical complete dataset at the uniform sampling interval Tu introduced in section 3a. In matrix form, the minimum-norm least squares solution for the Doppler spectrum estimate is given by Y = D, where

 
formula

D is the complex input time series vector, is the band-limiting diagonal matrix filter, is the discrete Fourier transform (DFT) matrix with elements F(k, l) = exp[−i2π(k − 1)(l − 1)/N], N is the length of the uniformly sampled (with spacing Tu) time series, and is a diagonal sampling matrix where only the diagonal indices corresponding to the presence of input data (in units of Tu) have unity elements, with all others zero. The superscript “T” denotes transpose, superscript “†” denotes conjugate transpose, and superscript “pinv” denotes pseudoinverse. This equation is identical to Eq. (6) in Chornoboy and Weber (1994), except that we opted to take the pseudoinverse instead of the strict inverse for better robustness.

To make the formulation specific for WSP, the elements of are defined as S[m, Q(m)] = 1 for k = 1, 2, …, M; otherwise S = 0. Here M is the number of time series data points available after the clutter filter (25, because the first and last points are dropped due to filter transients), Q(m) = (1, 10, 19, 28, 37, 46, 55, 64, 71, 78, 85, 92, 99, 106, 113, 120, 127, 134, 143, 152, 161, 170, 179, 188, 197), and N = 197.

At this point we diverge from Chornoboy and Weber (1994). Instead of iteratively adjusting the band-limiting filter until it converges on a dealiased solution, we employ a fixed and use the dealiased normalized mean frequency from section 3b to phase shift the input time series so that the first spectral moment is centered in the window defined by . This novel approach allows us to significantly reduce execution time by 1) precomputing offline and 2) eliminating the iterative dealiasing process. Kalogiros (2012) also uses frequency shifting in this context, but it is applied within an iterative solution like Chornoboy and Weber (1994).

We define filter to be 21 points in width: L(n, n) = 1 for n = 1, 2, …, 11 and 188, 189, …, N; otherwise L = 0. This corresponds to a spectral span of ± (about ±28 m s−1). This is wide enough to fit most weather spectra except those of some tornadoes. To center the mean Doppler shift of the input time series in this window, we multiply the elements of such that

 
formula

This phase-shifted vector is then multiplied by to arrive at the zero-Doppler-centered, band-limited Doppler spectrum vector, YShift = DShift. Next, the spectrum is shifted back to the correct location by using the normalized frequency s:

 
formula

where nShift = ROUND(sN), and “ROUND” denotes rounding to the nearest integer. The power spectrum vector is then computed as P = |Y|2/N.

Finally, because velocity dealiasing is limited to a realistic span, it is not necessary to keep all the spectral bins in storage. Only the number of central bins corresponding to the maximum dealiased velocity span plus the length of filter (corresponding to the nonzero diagonal elements) needs to be kept. We call this value NSpec.

d. Spectral buffer

Even with the ECPI, the ASR-9’s fast antenna rotation rate limits the number of pulses averaged per dwell to a small number compared to weather radars. Further averaging of the data is thus necessary for acceptable estimate error reduction. As explained in section 2, the legacy algorithm uses a recursive smoothing filter on lags 0 and 1 across scans for this purpose. A problem with this filter is that recovery from an anomalous value takes a long time. A median filter across scans will smooth the data just as well, but it will recover from an unrepresentative value more quickly. Although a median filter requires NMed times more memory, where NMed is the length of the median filter, the new DSP hardware easily allows for such an expansion. The spectral buffer of depth NMed is maintained for each beam (low and high) at every azimuth and range position. Currently, we have set NMed = 3, based on keeping the filter output variance similar to the variance for the legacy recursive filter with weight A = 0.5. It has been shown that the variance of power estimates from a square-law transfer function is almost the same for a uniform average of NFilt samples as for a first-order recursive filter if NFilt = (2 − A)/A (Zrnić 1977).

The probability distribution of weather I&Q signal power is expected to be exponential (Doviak and Zrnić 1993). Thus, the median of power over time will be different from the mean. As the filter length increases, the ratio of median to mean decreases from unity to an asymptote at ln(2). For NMed = 3, we calculated a normalization ratio of 0.83, which we use to adjust the median power spectra values to match the mean.

e. Range–Doppler filter

Moving clutter targets, such as aircraft and birds, can pull the computed spectrum to the wrong (i.e., not centered on the weather velocity) Doppler location. We have therefore adapted part of the moving clutter spectral filter (MCSF) that was developed for TDWR (Cho 2009). This operates like a point-target filter in the 2D range–Doppler domain.

Consider a range–Doppler array such that range is the vertical axis and Doppler is the horizontal axis. The filter runs along each Doppler spectral bin versus range, and the number of neighbors considered in the spectral dimension expands with range away from the range gate of interest. This expansion prevents filtering of cases where the spectral-range signature has a shape like “\,” “/,” “>,” or “<”, which are possible manifestations of wind shear. So, only features shaped like “—” or “•” that are indicative of unwanted nonmeteorological targets are filtered. To be more specific, at range gate j and spectral index n, compare the spectral power P(j, n) to PPeak(j − g, n) and PPeak(j + g, n), and mark P(j, n) for removal if it is stronger than both PPeaks by a threshold H(g). The gate difference g goes from 1 to 3, and the threshold values used currently are H(1) = 8 dB, H(2) = 10 dB, and H(3) = 17 dB. Removal is accomplished by linear interpolation over range. Term PPeak is the peak value of spectral power taken over bins kh to k + h (circularly wrapped around to the other side of the spectrum if needed). As stated earlier, this spectral neighborhood expands with range such that h = 1 for g = 1, h = 2 for g = 2, and h = 3 for g = 3.

The processing steps described in sections 3ae constitute an innovative algorithmic suite for generating clean Doppler spectral data for weather information extraction. To distinguish this collection of algorithms from the dual-beam velocity estimation procedure that follows, and for ease of reference, we dub it Moving Clutter Spectral Processing for Uneven-Sampled Data with Dealiasing (MCSPUDD).

f. Dual-beam velocity

The averaged and “cleaned” spectra are now used to generate reflectivity and velocity estimates. This process is straightforward for the single-beam quantities—they are computed in the usual way via R0 and R1 (e.g., Doviak and Zrnić 1993), where R0 is the total spectral power and

 
formula

where q(n) = n − floor(NSpec/2) − 1 and “FLOOR” denotes rounding down to the nearest integer. The velocity is then calculated as

 
formula

The dual-beam reflectivity is computed just as in the legacy algorithm, by taking the weighted average of the low- and high-beam reflectivities, where the weights depend on range. It is the estimation of dual-beam velocity where we apply a different technique.

As explained in the introduction, the fatness of the ASR-9 antenna beam in the vertical dimension is a detriment to measuring the velocity of a microburst outflow, which occurs as a thin layer over the ground surface. The legacy signal processing algorithm compensates for this deficiency by differencing the low- and high-beam autocorrelation lag-1 measurements [(1)] in an effort to focus the vertical resolution on the lower part of the low beam. This method, however, relies on an underdetermined weighting scheme.

A similar differencing operation can be performed in the Doppler spectral domain that does not involve weighting. It is also compatible with velocity dealiasing if the low- and high-beam spectra are dealiased prior to differencing as we are doing. Furthermore, it has been shown through analysis of real data that spectral differencing generates velocity estimates that are closer to measurements made of microburst outflows by a pencil-beam radar than estimates made using the autocorrelation differencing technique (Weber 1989). The spectral differencing algorithm was not initially implemented in WSP, because the processing and memory load would have been too much for the legacy DSP system. With the hardware upgrade we can now implement it.

The spectral differencing algorithm was developed independently by Weber and Noyes (1988) and Atlas (1989). Our implementation of it is straightforward. First, normalize the low- and high-beam power spectra so that their powers are equal. Second, subtract the normalized high-beam power spectrum from the normalized low-beam power spectrum; if the result is negative at a given spectral index, then set it to zero. If, anomalously, the resulting total power is zero, then replace the difference power spectrum with the low-beam power spectrum. Finally, the dual-beam velocity is computed from the difference power spectrum through (11) and (12).

The spectral differencing algorithm relies on the premise that in the presence of a vertical gradient in the horizontal velocity as in a microburst outflow, the separation in Doppler spectra corresponds to the separation in vertical space. Therefore, by subtracting the high-beam spectrum from the low-beam spectrum, the remaining Doppler spectrum more closely represents the velocity sensed by the bottom side of the low beam. As with any differencing operation, the price to pay is an increase in noise. The legacy algorithm compensates by applying a 3 × 3 (range × azimuth) median filter to the dual-beam velocity field in the postprocessing stage.

4. Algorithm performance evaluation

To test the performance of the new signal processing algorithms, we collected I&Q data with the FAA Academy’s ASR-9 WSP in Norman, Oklahoma. Data on three storms were captured on 12 and 25 June 2014. This ASR-9 is located at 35.4014°N, 97.6252°W, with an antenna altitude of 409 m. It operates at frequencies of 2720 and 2730 MHz—the latter channel was used during our data collection period. The mean values for the short- and long-PRI sets were 766 and 981 μs, respectively, corresponding to Nyquist velocities of ±36 and ±28 m s−1, respectively.

In the following data comparisons, the postprocessing step was not applied in order to exclude the effects of smoothing and censoring.

a. MCSPUDD

To provide an example comparison of the legacy signal processing with MCSPUDD, and at the same time illustrate how the new algorithm suite works, we focus on one radial of data. Figure 3 shows low-beam reflectivity data output from legacy processing with the plot zoomed into close range. Note the speckles in the reflectivity field that are likely due to unfiltered clutter, both moving and stationary. The white line indicates the particular azimuth (304°) radial that we will analyze. Figure 4 displays the corresponding radial velocity data from the legacy algorithm. There are patches of aliasing (red pixels against the blue background), some of which are traversed by our radial. The background flow is generally northwesterly.

Fig. 3.

Low-beam reflectivity estimates generated by the legacy WSP signal processing algorithm. The scan was conducted at 0657 UTC 12 Jun 2014 by the FAA Academy ASR-9. The white line indicates the azimuth position (304°) that will be discussed further.

Fig. 3.

Low-beam reflectivity estimates generated by the legacy WSP signal processing algorithm. The scan was conducted at 0657 UTC 12 Jun 2014 by the FAA Academy ASR-9. The white line indicates the azimuth position (304°) that will be discussed further.

Fig. 4.

As in Fig. 3, except the data shown are low-beam radial velocity.

Fig. 4.

As in Fig. 3, except the data shown are low-beam radial velocity.

Figure 5 shows band-limited Doppler spectral reconstruction versus range at azimuth 304°. Here, the spectral locations are based on undealiased velocity estimates. Thus, aliasing can clearly be seen at 6 and 7.3 km in range, corresponding to the red pixels seen in Fig. 4. The white line corresponds to mean velocities computed from the spectra using Eqs. (11) and (12).

Fig. 5.

Doppler spectral power vs range at azimuth 304° for the scan shown in Figs. 3 and 4. The band-limited spectral reconstructions were located within the larger spectral domain based on undealiased velocity estimates. The white line corresponds to mean velocities computed from the spectra.

Fig. 5.

Doppler spectral power vs range at azimuth 304° for the scan shown in Figs. 3 and 4. The band-limited spectral reconstructions were located within the larger spectral domain based on undealiased velocity estimates. The white line corresponds to mean velocities computed from the spectra.

In Fig. 6, the spectra were situated within the larger window based on the dealiased velocities (with false dealias correction), which rectified the aliased anomalies at 6 and 7.3 km.

Fig. 6.

As in Fig. 5, except the band-limited spectra were placed within the larger spectral domain based on dealiased velocity estimates.

Fig. 6.

As in Fig. 5, except the band-limited spectra were placed within the larger spectral domain based on dealiased velocity estimates.

In Fig. 7, the three-scan median filter was applied, which provided some smoothing. The spectral scatter around 5 km in range was due to variability over time; thus, the temporal median filter was able to generate a more stable spectrum in those range gates.

Fig. 7.

As in Fig. 6, except the spectra were additionally processed with a three-point median filter across consecutive scans.

Fig. 7.

As in Fig. 6, except the spectra were additionally processed with a three-point median filter across consecutive scans.

Finally, the 2D range–Doppler moving clutter filter was added in Fig. 8, successfully removing the Doppler anomalies at ranges 0.6 and 2.5 km, resulting in very clean spectral and mean velocity profiles. Note that each processing step contributed to the improvement of data quality in this example.

Fig. 8.

As in Fig. 7, except that the moving clutter range–Doppler filter was applied.

Fig. 8.

As in Fig. 7, except that the moving clutter range–Doppler filter was applied.

Figures 9 and 10 show the MCSPUDD output for the same scan that was used in the legacy processing output of Figs. 3 and 4. Much of the speckle in the reflectivity has disappeared, and the velocity field is not as marred with aliasing and inconsistent patches. Velocity irregularity in the top part of Fig. 10 corresponding to areas of low reflectivity (SNR) is to be expected. This is not a cherry-picked example—it is quite representative of most scans. Successful velocity dealiasing was seen consistently from scan to scan and during other periods, and moving clutter targets were effaced effectively.

Fig. 9.

As in Fig. 3, except data were processed using the new signal processing algorithm (MCSPUDD).

Fig. 9.

As in Fig. 3, except data were processed using the new signal processing algorithm (MCSPUDD).

Fig. 10.

As in Fig. 4, except data were processed using the new signal processing algorithm (MCSPUDD).

Fig. 10.

As in Fig. 4, except data were processed using the new signal processing algorithm (MCSPUDD).

b. Spectral differencing dual-beam velocity

To compare the dual-beam velocity estimation performance of the old and new algorithms, we needed a source of “truth.” Fortunately, the FAA Academy’s ASR-9 is located very close to two TDWRs—the Academy (OEX) and Program Support Facility (PSF) units. The OEX TDWR is at 35.4053°N, 97.6252°W, 0.64 km on azimuth 312° from the ASR-9, with an antenna altitude of 416 m. The PSF TDWR is at 35.3927°N, 97.6287°W, 1.25 km on azimuth 219° from the ASR-9, with an antenna altitude of 424 m. The OEX and PSF TDWRs are available for experiments on an ad hoc basis when they are not being used for instructional or testing purposes. They are both close enough in location and altitude to the ASR-9 that direct comparison of data is possible, simply by shifting the resampled data grids. TDWR operates at C band and has a conical antenna beam of width 0.55°.

At the default installation tilt (0°) of the ASR-9 antenna, its low-beam elevation pattern reaches its maximum at 2°. The academy ASR-9 antenna is installed with a tilt of 1.5°, so the low-beam nose points at 3.5°, with the high-beam maximum at 7.5° (Fig. 11). For reference, the hazardous volume scan elevation angles for the OEX TDWR are also shown in Fig. 11. The 12 June storm data were collected by the OEX TDWR in this mode. Thus, the WSP dual-beam velocity should be compared to the 2.6° elevation angle data from the OEX TDWR. The 25 June data were collected by the PSF TDWR running in monitor volume scan mode for which the low-elevation angles were 0.3°, 1°, 2.5°, 6.1°, 11°, and 15.9°. So, the 2.5° elevation angle data from the PSF TDWR should be used for dual-beam velocity comparison, as it lies on the bottom side of the ASR-9 low-beam pattern. The antenna azimuth rotation rates were 26° s−1 for the OEX hazard scan at 2.6° elevation and 19° s−1 for the PSF monitor scan at 2.5° elevation.

Fig. 11.

Academy ASR-9 two-way antenna patterns vs elevation angle. Indicated on the right margin are the academy TDWR hazardous volume scan elevation angles.

Fig. 11.

Academy ASR-9 two-way antenna patterns vs elevation angle. Indicated on the right margin are the academy TDWR hazardous volume scan elevation angles.

Theoretically, the weighted average of TDWR data from multiple elevation cuts would provide a better comparison to data taken with the ASR-9’s fan beam. In practice, however, we have found this approach to be problematic. First, because the TDWR data at each elevation angle has invalid data flags at different locations, combining data from multiple elevation scans dramatically decreases the amount of available valid data. Second, because the TDWR’s hazardous volume scan interleaves scans nonmonotonically in elevation, the cuts that are averaged are more spread out in time. Hence, we decided to make the comparisons at single elevation scans of the TDWR.

To prepare for comparison, we resampled the range–azimuth velocity data from all radars in Cartesian coordinates at 100-m resolution. We then shifted the TDWR data—500 m eastward and 400 m southward for OEX, 800 m eastward and 1000 m northward for PSF. We restricted the comparison domain to the central 16 km × 16 km square, since the WSP microburst detection algorithm, for which the dual-beam velocity is used, only extends to 16-km range. Pixels where the TDWR data were flagged invalid were eliminated. The comparisons were made at the ASR-9 scan times that were closest to the TDWR scan times.

For a visual comparison, Figs. 1214 show the TDWR radial velocity, the legacy WSP dual-beam velocity, and the new spectral-differencing WSP dual-beam velocity. The inbound magnitudes in the northwest sector (the blues) are better matched by the new algorithm, as are the outbound magnitudes (yellows and greens) to the southwest. The contrast in performance can be observed more explicitly by taking the absolute value of the velocity differences (Figs. 15 and 16). It can be seen readily that overall the new algorithm did a better job of matching the TDWR data.

Fig. 12.

Radial velocity plot generated by the academy TDWR at 1008:08 UTC 12 Jun 2014. The elevation angle was 2.6°. Censored areas due to range aliasing, velocity dealias failure, and low signal power relative to noise or clutter are shown in white.

Fig. 12.

Radial velocity plot generated by the academy TDWR at 1008:08 UTC 12 Jun 2014. The elevation angle was 2.6°. Censored areas due to range aliasing, velocity dealias failure, and low signal power relative to noise or clutter are shown in white.

Fig. 13.

Dual-beam velocity plot generated for the academy ASR-9 WSP data at 1008:11 UTC 12 Jun 2014 with the legacy signal processing algorithm.

Fig. 13.

Dual-beam velocity plot generated for the academy ASR-9 WSP data at 1008:11 UTC 12 Jun 2014 with the legacy signal processing algorithm.

Fig. 14.

As in Fig. 13, except with the data processed using the new signal processing algorithm.

Fig. 14.

As in Fig. 13, except with the data processed using the new signal processing algorithm.

Fig. 15.

Plot of the absolute value of the difference between the TDWR velocity data of Fig. 12 and the legacy WSP dual-beam velocity data of Fig. 13.

Fig. 15.

Plot of the absolute value of the difference between the TDWR velocity data of Fig. 12 and the legacy WSP dual-beam velocity data of Fig. 13.

Fig. 16.

Plot of the absolute value of the difference between the TDWR velocity data of Fig. 12 and the new spectral-differencing WSP dual-beam velocity data of Fig. 14.

Fig. 16.

Plot of the absolute value of the difference between the TDWR velocity data of Fig. 12 and the new spectral-differencing WSP dual-beam velocity data of Fig. 14.

Using the root-mean-square (RMS) error of the difference with respect to the TDWR velocity data as the performance metric (averaged over all pixels in the comparison domain), we plot the comparison in Fig. 17. The new dual-beam velocity algorithm consistently outperformed the legacy algorithm through three storms. The RMS errors averaged over all events were 8.2 m s−1 for the legacy algorithm and 7.0 m s−1 for the new algorithm, which is a 15% reduction in estimation error.

Fig. 17.

RMS error relative to TDWR velocity data for the legacy and new WSP dual-beam velocity algorithms. Each case number represents a scan comparison averaged over all valid pixels in the central 16 km × 16 km domain. The TDWR elevation angle was 2.6° during the 12 Jun cases and 2.5° during the 25 Jun cases.

Fig. 17.

RMS error relative to TDWR velocity data for the legacy and new WSP dual-beam velocity algorithms. Each case number represents a scan comparison averaged over all valid pixels in the central 16 km × 16 km domain. The TDWR elevation angle was 2.6° during the 12 Jun cases and 2.5° during the 25 Jun cases.

To show that the computed WSP low- and high-beam velocities were, indeed, closest to the velocities from the expected TDWR elevation angle scans, we plot the overall RMS errors versus elevation angle in Fig. 18. For the low beam, the error was smallest relative to the 2.6° OEX elevation angle (2.5° for PSF), and for the high beam the error was smallest at the 5.3° OEX elevation angle (6.1° for PSF), which is consistent with the antenna elevation beam plot of Fig. 11. The new processing algorithm outperformed the legacy algorithm for these beams also (8% and 15% error reductions for the low and high beams, respectively), which we attribute to the cleaning and dealiasing action of MCSPUDD. Note that by comparing the velocity error profiles of the low-beam and dual-beam cases in Fig. 18, we can discern the “center of gravity” of the effective observed elevation angle for the latter being lower than the former; that is, the low-beam error was greater at 1° than at 5.3°, whereas the opposite was true for the dual-beam error. If the TDWR scan angles were more densely distributed, then the effect should have been even clearer.

Fig. 18.

RMS error relative to TDWR velocity data for the legacy and new WSP velocity algorithms computed over all cases for the (left) low beam, (middle) high beam, and (right) dual beam. The elevation angles shown are for the OEX TDWR; the 25 Jun data from the PSF TDWR were taken at slightly different angles for the top two positions shown (2.5° and 6.1°).

Fig. 18.

RMS error relative to TDWR velocity data for the legacy and new WSP velocity algorithms computed over all cases for the (left) low beam, (middle) high beam, and (right) dual beam. The elevation angles shown are for the OEX TDWR; the 25 Jun data from the PSF TDWR were taken at slightly different angles for the top two positions shown (2.5° and 6.1°).

Note that the dual-beam velocity error was slightly larger than the low-beam error at 2.6°. This is due to the dual-beam velocity estimate resulting from a differencing operation, which amplifies the noise. As alluded to before, this increase in variance is alleviated in the postprocessing step by applying a spatial median filter, which was not included in our analysis.

Also note that reflectivity estimate comparisons with TDWR data were not attempted, because the TDWR, operating at C band, suffers significantly more attenuation through storms as well as more range-aliased contamination from distant weather than does the ASR-9. Therefore, the TDWR is not a reliable standard for reflectivity estimate comparison.

5. Summary discussion

The hardware upgrade of the ASR-9 WSP DSP will allow an enhancement of the signal processing algorithms. This paper showed that a suite of new algorithms dubbed MCSPUDD can improve reflectivity and radial velocity estimates by filtering out isolated moving clutter targets and dealiasing velocity. Comparison with TDWR data showed that low- and high-beam velocity estimation accuracy are improved by 8% and 15%, respectively. Furthermore, replacement of the autocorrelation-lag-1-based dual-beam velocity calculation by a spectral-differencing approach can extend velocity dealiasing to the dual-beam product, and it can increase estimate accuracy. Comparison with TDWR data showed that dual-beam velocity accuracy is improved by 15%.

The moving clutter filtering aspect of MCSPUDD is especially beneficial for hazardous terminal weather surveillance, since there are usually many aircraft in the same airspace that is being covered for wind shear alerts. Because of its vertical fan beam, the ASR-9 is also more likely to have an aircraft (or birds and bats) in its beam at any given range–azimuth cell than a TDWR with its narrow pencil beam.

Although the TDWRs were close enough to the ASR-9 for excellent geometric coincidence in the horizontal dimension, the rather coarse sampling in elevation angle by the TDWR volume scans (Fig. 11) made for imperfect data comparisons. The dual-beam velocity comparison had a good corresponding TDWR scan angle (on the bottom side of the ASR-9 low-beam pattern), but the low- and high-beam comparisons did not. Ideally, the TDWR elevation cuts would have matched the ASR-9 elevation beam pattern maxima. Therefore, the 8% and 15% velocity error reductions for the low and high beams, respectively, are likely to be underestimates. For future algorithm performance evaluations, it would be beneficial to design custom TDWR volume scans that sample only low-elevation angles with closer spacing. For reflectivity estimate comparisons, an S-band pencil-beam radar—that is, a NEXRAD—located nearby would be needed in order to eliminate the attenuation and range-aliasing problems associated with the C-band TDWR.

In the long term, a second weather receiver channel added to the WSP would enable simultaneous reception of high- and low-beam I&Q data, which would allow coherent processing techniques like interferometry. This could potentially provide a better way to measure near-surface winds with the ASR-9, but at this time there are no plans to add the extra hardware due to cost considerations.

Acknowledgments

This work was sponsored by the Federal Aviation Administration under U.S. Air Force Contract FA8721-05-C-0002. The opinions, interpretations, conclusions, and recommendations are those of the author and are not necessarily endorsed by the U.S. government. I would like to thank Cory West and Alok Gautam of the FAA Program Support Facility (PSF) in Oklahoma City, Oklahoma, for collecting the ASR-9 WSP I&Q data and Pete Smith (also PSF) for obtaining the TDWR base data used in this research; John Morgan for implementing the algorithms in real-time code and making changes to reduce execution time; and Geoff Stilwell for data transfer and archival assistance. I would also like to acknowledge the FAA program manager, Steve Kim, for supporting this work.

REFERENCES

REFERENCES
Atlas
,
D.
,
1989
: The detection of low level windshear with airport surveillance radar. Preprints, Third Int. Conf. on the Aviation Weather System, Anaheim, CA, Amer. Meteor. Soc.,
21
24
.
Cho
,
J. Y. N.
,
2005
:
Multi-PRI signal processing for the Terminal Doppler Weather Radar. Part II: Range–velocity ambiguity mitigation
.
J. Atmos. Oceanic Technol.
,
22
,
1507
1519
, doi:.
Cho
,
J. Y. N.
,
2009
: Moving clutter spectral filter for Terminal Doppler Weather Radar. 34th Conf. on Radar Meteorology, Williamsburg, VA, Amer. Meteor. Soc., P5.2. [Available online at https://ams.confex.com/ams/34Radar/techprogram/paper_155381.htm.]
Cho
,
J. Y. N.
, and
E. S.
Chornoboy
,
2005
:
Multi-PRI signal processing for the Terminal Doppler Weather Radar. Part I: Clutter filtering
.
J. Atmos. Oceanic Technol.
,
22
,
575
582
, doi:.
Chornoboy
,
E. S.
,
1993
: Clutter filter design for multiple-PRT signals. Preprints, 26th Conf. on Radar Meteorology, Norman, OK, Amer. Meteor. Soc.,
235
237
. [Available online at http://www.ll.mit.edu/mission/aviation/publications/publication-files/ms-papers/Chornoboy_1993_RM_MS-10254_WW-18698.pdf.]
Chornoboy
,
E. S.
, and
M. E.
Weber
,
1994
: Variable-PRI processing for meteorologic Doppler radars. Record of the 1994 IEEE National Radar Conference, IEEE,
85
90
. [Available online at http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=328103.]
Ding
,
C.
,
D.
Pei
, and
A.
Salomaa
,
1996
:
Chinese Remainder Theorem: Applications in Computing, Coding, Cryptography.
World Scientific Publishing, 213 pp.
Doviak
,
R. J.
, and
D. S.
Zrnić
,
1993
:
Doppler Radar and Weather Observations.
Academic Press, 562 pp.
Kalogiros
,
J.
,
2012
:
Least squares reconstruction of Doppler radar spectra for irregular PRT
.
J. Atmos. Oceanic Technol.
,
29
,
1744
1756
, doi:.
Michelson
,
M.
,
W. W.
Shrader
, and
J. G.
Wieler
,
1990
:
Terminal Doppler Weather Radar
.
Microwave J.
,
33
,
139
148
.
Taylor
,
J. W.
, Jr.
, and
G.
Brunins
,
1985
:
Design of a new airport surveillance radar (ASR-9)
.
Proc. IEEE
,
73
,
284
289
, doi:.
Torres
,
S. M.
,
Y. F.
Dubel
, and
D. S.
Zrnić
,
2004
:
Design, implementation, and demonstration of a staggered PRT algorithm for the WSR-88D
.
J. Atmos. Oceanic Technol.
,
21
,
1389
1399
, doi:.
Trunk
,
G.
, and
S.
Brockett
,
1993
: Range and velocity ambiguity reduction. Record of the 1993 IEEE National Radar Conference, IEEE Aerospace and Electronic Systems Society,
146
149
.
Weber
,
M. E.
,
1989
: Dual-beam autocorrelation based wind estimates from airport surveillance radar signals. MIT Lincoln Laboratory Project Rep. ATC-167, DOT Rep. DOT/FAA/PS-89/5, 60 pp. [Available online at http://www.ll.mit.edu/mission/aviation/publications/publication-files/atc-reports/Weber_1989_ATC-167_WW-15318.pdf.]
Weber
,
M. E.
,
2002
: ASR-9 weather systems processor (WSP) signal processing algorithms. MIT Lincoln Laboratory Project Rep. ATC-255, 53 pp. [Available online at http://www.ll.mit.edu/mission/aviation/publications/publication-files/atc-reports/Weber_2002_ATC-255_WW-15318.pdf.]
Weber
,
M. E.
, and
T. A.
Noyes
,
1988
: Low-altitude wind shear detection with airport surveillance radars: Evaluation of 1987 field measurements. MIT Lincoln Laboratory Project Rep. ATC-159, DOT Rep. DOT/FAA/PS-88/10, 120 pp. [Available online at http://www.ll.mit.edu/mission/aviation/publications/publication-files/atc-reports/Weber_1988_ATC-159_WW-15318.pdf.]
Weber
,
M. E.
, and
M. L.
Stone
,
1995
:
Low altitude wind shear detection using airport surveillance radars
.
IEEE Aerosp. Electron. Syst. Mag.
,
10
,
3
9
, doi:.
Weber
,
M. E.
,
J. A.
Cullen
,
S. W.
Troxel
, and
C. A.
Meuse
,
1996
: ASR-9 weather system processor (WSP): Wind shear algorithms performance assessment. MIT Lincoln Laboratory Project Rep. ATC-247, 42 pp. [Available online at http://www.ll.mit.edu/mission/aviation/publications/publication-files/atc-reports/Weber_1996_ATC-247_WW-15318.pdf.]
Wilson
,
F. W.
, and
R. H.
Gramzow
,
1991
: The redesigned low level wind shear alert system. Preprints, Fourth Int. Conf. on Aviation Weather Systems, Paris, France, Amer. Meteor. Soc., 24–26.
Zrnić
,
D. S.
,
1977
:
Mean power estimation with a recursive filter
.
IEEE Trans. Aerosp. Electron. Syst.
,
AES-13
,
281
289
, doi:.